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Summary 

Work  was  conducted  on  two  projects  related  to  real-time  control  of  shear  flows.  In  the 
first,  two-dimensional  unsteady  simulations  of  the  development  of  the  wake  behind  a  circular 
cylinder  impulsively  started  into  rotatory  and  rectilinear  motion  were  performed.  This 
simulation  code  is  now  serving  as  a  testbed  for  the  development  of  open-  and  closed-loop 
strategies  for  control  of  the  lift/drag  ratio  as  well  as  the  suppression  of  vortex  shedding.  The 
code  has  been  checked  by  comparison  to  earlier  computational  and  experimental  work.  In  the 
second  project,  a  technique  was  developed  to  extract  both  velocity  components  in  a  two- 
dimensional  incompressible  flow  from  measurements  of  a  single  scalar  (temperature  or 
concentration),  and  all  three  velocity  components  in  a  three-dimensional  incompressible  flow 
from  measurements  of  two  scalars.  The  technique  is  applicable  to  steady  or  unsteady,  laminar 
or  turbulent  flows.  A  key  advantage  over  particle  image  velodmetry  and  other  multi-point 
techniques  is  that  our  method  uses  full-field  optical  measurements,  so  that  spatial  resolution 
is  not  limited  by  particle  size  and  loading  restrictions. 

1 .  Development  of  the  Wake  behind  a  Circular  Cylinder  Impulsively 
Started  Into  Rotatory  and  Translational  Motion 

We  (Chen,  Ou  &  Pearlstein  1993)  have  completed  a  computational  study  of  the  temporal 
development  of  two-dimensional  incompressible  flow  generated  by  a  circular  cylinder 
impulsively  started  into  steady  rotatory  and  rectilinear  motion  at  Re  -  200  (based  on  the 
cylinder  diameter  2a  and  the  magnitude  U  of  the  rectilinear  velocity).  We  have  used  an  explicit 
finite-difference/pseudo-spectral  technique  and  a  new  implementation  of  the  Biot-Savart  law 
to  integrate  a  velocity/vorticity  formulation  of  the  Navier-Stokes  equations.  Results  are 
presented  for  the  four  angular/rectilinear  speed  ratios  a-  fla/U  (where  W  is  the  angular 
speed)  considered  experimentally  by  Coutanceau  &  M6nard  (1985).  For  a  £  1,  extension  of 
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the  computations  to  larger  dimensionless  times  than  possible  in  the  work  of  Coutanceau  & 
M6nard  or  considered  computationally  by  Badr  &  Dennis  (1985)  allows  for  a  more  complete 
discussion  of  the  long-time  development  of  the  wake.  We  discuss  several  aspects  of  the  vortex 
kinematics  and  dynamics  not  revealed  by  the  earlier  work,  in  which  vortex  cores  were 
identified  from  frame-dependent  streamline  and  streamfunction  information,  rather  than  from 
the  frame  invariant  vorticity  distribution,  as  in  the  present  work.  Our  results  indicate  that  at 
Re  -  200,  vortex  shedding  does  indeed  occur  for  a  -  3.25  (and  possibly  for  larger  a),  in 
contrast  to  the  conclusion  of  Coutanceau  &  Mdnard.  Moreover,  the  shedding  process  is  very 
different  from  that  associated  with  the  usual  Kdrmdn  vortex  street  for  a  «  0.  Specifically, 
consecutive  vortices  can  be  shed  from  one  side  of  the  cylinder  and  be  of  the  same  sense,  in 
contrast  to  the  nonrotating  case,  in  which  mirror-image  vortices  of  opposite  sense  are  shed 
alternately  from  opposite  sides  of  the  cylinder. 

The  code  is  now  being  used  in  the  development  of  open-  and  dosed  loop  strategies  for 
controlling  the  lift/drag  ratio  (Ou  &  Bums)  and  suppressing  vortex  shedding  (Tsai,  Bentsman 
&  Pearistein).  A  copy  of  a  paper  published  in  J.  Fluid  Mech.  is  included  as  Appendix  A,  and  a 
reprint  is  endosed. 

2.  Determination  of  Two-  and  Three-Dimensional  Incompressible  Velocity 
Fields  from  Measurements  of  Passive  or  Reactive  Scalars 

We  have  developed  a  technique  for  determining  the  velocity  fields  in  two-  and  three- 
dimensional  incompressible  flows  from  full-field  measurements  of  one  or  two  passive  or 
reactive  scalars.  Our  technique  is  based  on  the  recognition  that  if  the  transport  equation(s)  for 
the  scaiar(s)  is  (are)  linear,  then  it  (they)  and  the  incompressible  continuity  equation 

V.u=0,  (1) 


constitute  two  (three)  linear  equations  in  two  (three)  unknown  velocity  components.  The 
method  is  fully  described  in  Appendix  B,  which  is  a  preprint  of  a  manuscript  submitted  to 
Physics  of  Fluids. 

in  n  dimensions,  we  reduce  (1)  and  n-  scalar  transport  equations  to  n  uncoupled  linear 
hyperbolic  equations  in  the  n  velocity  components.  Given  the  scalar  field  data,  solution  of  this 
systems  is  made  unique  by  spedfication  of  appropriate  boundary  conditions.  Boundary  values 
of  the  unknown  velocity  components,  rather  than  the  scalars  are  required.  We  note  that  only 
spatial  derivatives  of  the  velocity  components  appear  in  these  equations.  Thus  no  initial  data  on 
the  velocity  field  are  required.  A  desirable  consequence  of  the  fact  that  the  hyperbolic 
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equations  are  not  initial  value  problems  is  that  errors  in  the  determination  of  the  velocity  field 
do  not  grow  in  time. 

In  two  dimensions,  the  equations  are  singular  at  any  point  at  which  the  gradient  of  the 
scalar  vanishes,  and  in  three  dimensions  they  are  singular  at  points  for  which  one  or  more 
components  of  the  cross  product  of  the  gradients  of  the  scalars  vanish.  We  have  shown  that 
these  singularities  are  completely  removable  (in  the  sense  that  the  extracted  velocity  field  is 
smooth  and  can  be  uniquely  determined)  unless  the  gradient  of  the  scalar  vanishes  identically 
over  an  area  (in  two  dimensions)  or  the  cross  product  of  the  gradients  of  the  scalars  vanishes 
identically  over  a  volume  (in  three  dimensions).  In  that  case,  the  scalar  field  does  not  provide 
enough  information  to  determine  the  velocity  field,  and  any  method  necessarily  fails. 

We  have  also  analyzed  three  methods  which  claim  to  extract  the  velocity  field  (for 
incompressible  or  compressible  flows)  from  measurements  of  a  single  scalar.  For  one,  the 
"iterative  inversion"  technique  (Dahm  et  al.  1992),  it  is  shown  that  if  convergence  does  occur, 
then  that  result  must  be  a  consequence  of  spatial  discretization,  since  in  the  limit  of  "full" 
resolution,  it  can  be  proved  that  the  proposed  iteration  cannot  converge.  For  the  second, 
"variational"  technique  of  Tokumaru  &  Dimotakis  (1994)  which  seeks  the  "smoothest" 
velocity  field,  we  point  out  that  although  the  correct  flow  is  not  always  the  smoothest  one 
consistent  with  the  transport  data,  there  are  likely  to  be  situations  in  which  the  correct 
velocity  field  and  the  smoothest  velocity  field  do  not  differ  by  very  much.  Finally,  we  show 
that,  for  a  second  "variational"  technique,  the  criteria  imposed  on  the  velocity  field  are  not 
satisfied  unless  either  the  flow  is  steady  and  inertial  terms  are  negligible,  or  the  viscous  terms 
can  be  written  as  the  gradient  of  a  scalar  potential.  Neither  of  these  situations  is  likely  to 
obtain  in  flows  of  interest. 

Finally,  we  have  completed  a  two-dimensional  proof-of-concept  study  (Carpenter  1993) 
showing  that  the  velocity  held  can  be  extracted  from  a  computed  scalar  in  a  diverging  channel 
flow.  The  extracted  velocity  field  is  nonparallel,  and  fully  two-dimensional.  The  results 
clearly  show  that  the  solenoidal  velocity  field  can  be  extracted  from  a  single  scalar  in  two 
dimensions.  A  detailed  study  of  the  effects  of  noise  and  spatial  resolution  is  underway. 

Joint  work  with  Professor  An  Glezer  is  underway  to  implement  this  method  in  the  plane 
mixing  layer  facility  at  Georgia  Tech,  with  an  eye  toward  real-time  sensing  and  control. 
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Development  of  the  wake  behind  a  circular 
cylinder  impulsively  started  into  rotatory  and 
rectilinear  motion 
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The  temporal  development  of  two-dimensional  viscous  incompressible  flow  generated 
by  a  circular  cylinder  impulsively  started  into  steady  rotatory  and  rectilinear  motion  at 
Re  =  200  (based  on  the  cylinder  diameter  2 a  and  the  magnitude  U  of  the  rectilinear 
velocity)  is  studied  computationally.  We  use  an  explicit  finite-difference/pseudospectral 
technique  and  a  new  implementation  of  the  Biot-Savart  law  to  integrate  a 
velocity/vorticity  formulation  of  the  Navier-Stokes  equations.  Results  are  presented 
for  the  four  angular: rectilinear  speed  ratios  a  -  Qa/U  (where  Q  is  the  angular  speed) 
considered  experimentally  by  Coutanceau  &  Menard  (1985).  For  at  ^  1,  extension  of 
the  computations  to  dimensionless  times  larger  than  achieved  either  in  the  experimental 
work  or  in  the  computations  of  Badr  &  Dennis  (1985)  allows  for  a  more  complete 
discussion  of  the  temporal  development  of  the  wake.  Using  the  frame-invariant 
vorticity  distribution,  we  discuss  several  aspects  of  the  vortex  kinematics  and  dynamics 
not  revealed  by  the  earlier  work,  in  which  vortex  cores  were  identified  from  frame- 
dependent  streamline  and  streamfunction  information.  Consideration  of  the  flow  in  the 
absence  of  sidewalls  confirms  the  artifactual  nature  of  the  trajectory  of  the  first  vortex 
reported  by  Coutanceau  &  Menard  for  a  =  3.25.  For  a  greater  than  unity  (the  largest 
value  considered  by  Badr  &  Dennis),  our  results  indicate  that  at  Re  =  200  shedding  of 
more  than  one  vortex  does  indeed  occur  for  at  =  3.25  (and  possibly  for  larger  a),  in 
contrast  to  the  conclusion  of  Coutanceau  &  Menard.  Moreover,  the  shedding  process 
is  very  different  from  that  associated  with  the  usual  Karman  vortex  street  for  at  =  0. 
Specifically,  consecutive  vortices  can  be  shed  from  one  side  of  the  cylinder  and  be  of 
the  same  sense,  in  contrast  to  the  non-rotating  case,  in  which  mirror-image  vortices  of 
opposite  sense  are  shed  alternately  from  opposite  sides  of  the  cylinder.  The  results  are 
discussed  in  relation  to  the  possibility  of  suppressing  vortex  shedding  by  open-  or 
closed-loop  control  of  the  rotation  rate. 


t  Present  Address:  General  Electric  Aircraft  Engines,  1  Neumann  Way,  Cincinnati,  OH  45215, 
USA. 

t  Present  Address:  Interdisciplinary  Center  for  Applied  Mathematics.  Virginia  Polytechnic 
Institute  and  State  University,  Blacksburg,  VA  24061,  USA. 
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1.  Introduction 

Flow  past  a  rotating  circular  cylinder  is  a  prototypical  problem  in  the  study  of 
unsteady  flow  separation  (Telionis  1981).  It  is  also  of  considerable  practical  importance 
in  boundary-layer  control  on  airfoils  (cf.  Tennant,  Johnson  &  Krothapalli  1976  and 
Modi,  Mokhtarian  &  Yokomizo  1990),  and  in  lift  enhancement  schemes  employing  the 
Magnus  effect  (Swanson  1961).  Rotation  of  all  or  part  of  a  body  may  also  have 
applications  in  active  or  feedback  control  of  vortex  shedding,  with  important 
consequences  for  wake  modification  and  the  reduction  of  flow-induced  vibration. 

In  this  work,  we  describe  the  development  of  the  two-dimensional  flow  generated  by 
a  circular  cylinder  of  radius  a  started  impulsively  into  combined  steady  rotatory  and 
rectilinear  motion,  with  angular  speed  Q  about  its  axis  and  rectilinear  speed  U  normal 
to  its  generators.  The  fluid  is  taken  to  be  at  rest  initially.  The  two  parameters  governing 
the  development  of  the  flow  are  the  Reynolds  number,  defined  by  Re  =  2aU/v,  where 
v  is  the  kinematic  viscosity,  and  the  ratio  of  rotatory  to  rectilinear  speeds,  defined  by 
a  =  Qa/U. 

Experimental  studies  of  the  nominally  two-dimensional  flow  past  a  circular  cylinder 
undergoing  steady  rotatory  and  rectilinear  motion  have  been  conducted  by  Prandtl 
(Prandtl  1925;  Prandtl  &  Tietjens  1934),  Taneda  (1977,  1980),  Koromilas  &  Telionis 
(1980),  Diaz  et  al.  (1983),  Werle  (1984),  and  Kimura,  Tsutahara  &  Wang  (1992).  The 
most  detailed  work  is  that  of  Coustanceau  &  Menard  (1985)  and  Badr  et  al.  (1990),  in 
which  papers  an  excellent  summary  of  earlier  work  can  be  found.  On  the  basis  of  their 
experiments  (primarily  at  Re  =  200,  but  including  results  for  Re  as  high  as  1000), 
Coutanceau  &  Menard  (1985)  concluded  that  a  (modified)  Karman  vortex  street 

disappears  completely  for  a.  greater  than  a  certain  limiting  value  a.L.  The  value  of  aL  has 
been  found  to  be  not  very  dependent  on  the  Reynolds  number  and  to  be  about  2.  For 
a  >  aL  no  other  eddy  is  created  after  £,  (the  first  eddy  formed)  during  the  time  of  the 
observations,  so  that  the  eddy  street  must  have  been  destroyed. 

They  found  this  conclusion  to  be  consistent  with  the  earlier  experimental  work  at 
higher  Re  by  Prandtl,  Diaz  et  al.,  and  Werle  and  proposed,  as  a  simple  physical 
explanation  for  the  disappearance  of  the  Karn’in  vortex  street  at  sufficiently  high 
values  of  at,  that 

for  low  values  of  a,  eddies  would  be  alternately  shed  on  each  side  of  the  cylinder  to  form 
a  Benard-Karman  street,  as  for  the  pure  translation  (a  =  0).  But  the  eddies  on  the  side 
moving  in  the  direction  of  the  rotation  decrease  progressively  when  at  increases  and  then 
disappear  completely.  Thus  it  was  found  that  the  Benard-Karman  structure  begins  to 
deteriorate  as  soon  as  the  peripheral  velocity  becomes  greater  than  the  free-stream 
velocity  (giving  rise  to  a  zigzag  oscillating  wake)  and  finally  disappears  for  a  >  2.5. 

When  one  examines  the  evidence  for  these  statements,  one  finds  that  it  is  not 
overwhelmingly  strong  at  lower  values  of  Re,  particularly  for  the  critical  value  of  a  and 
its  dependence  on  Re.  For  Re  =  9000  and  several  values  of  a,  Diaz  et  al.  (1983)  made 
hot-wire  measurements  of  the  streamwise  velocity,  and  computed  its  autocorrelation. 
They  found  that  for  a  =  0  and  1,  the  velocity  autocorrelations  were  very  similar, 
approximately  periodic,  and  had  local  maxima  separated  by  a  time  corresponding  to 
the  nominal  Strouhal  frequency.  For  a  =  1.5  and  2,  the  autocorrelation  function  was 
progressively  reduced.  Diaz  et  al.  (1983)  did  indeed  conclude  that  ‘for  peripheral 
velocities  up  to  the  value  of  the  free-stream  velocity,  a  distinct  Karman  vortex  activity 
exists  within  the  wake,  whereas  for  greater  peripheral  velocities,  the  Karman  activity 
deteriorates  and  disappears  for  values  in  excess  of  twice  the  free-stream  velocity.’  On 
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the  other  hand,  for  Re  »  3300  Werle  (1984)  noted  that  for  values  of  a  in  excess  of  that 
at  which  separation  is  eliminated,  ‘  when  the  tangential  velocity  increases  further,  the 
cylinder  finally  entrains  an  entire  layer  of  relatively  turbulent  fluid  in  its  rotation.  More 
or  less  periodic  instabilities  then  appear.'  From  this,  it  is  not  clear  whether  vortex 
shedding  is  really  suppressed  by  rotation  at  Re  %  3300.  At  Re  =  10s,  Badr  et  al.  (1990) 
found  experimentally  and  computationally  that  the  second  vortex  was  formed  much 
later  for  a  =  2  than  for  a  =  1,  and  was  also  much  weaker.  For  a  =  3.  their  experiments 
and  computations  showed  that  the  first  two  vortices  formed  were  of  the  same 
rotational  sense,  and  that  one  of  the  vortices  is  shed  downstream,  while  the  other  ‘  is 
washed  down  to  the  frontal  part  of  the  cylinder  and  disappears’.  For  larger 
dimensionless  times,  their  two-dimensional  computations  show  that  no  additional 
vortices  are  formed,  and  the  computed  flow  approaches  a  steady  state.  In  their 
experiments,  ‘three-dimensional  and  instability  effects  become  more  pronounced, 
especially  in  the  wake’.  At  the  lower  Re  in  the  experiments  of  Coutanceau  &  Menard 
(1985),  the  towing  tank  used  allowed  observations  to  be  made  over  only  a  very  limited 
range  of  dimensionless  time. 

The  issues  of  whether  the  Karman  vortex  street  is  destroyed  and  vortex  shedding  is 
suppressed  are  of  considerable  practical  interest  from  the  standpoint  of  wake 
modification  and  the  reduction  of  flow-induced  vibration.  In  particular,  it  is  of  interest 
to  determine  whether,  for  a  given  Re ,  there  is  a  value  of  aL  beyond  which  (two-  or 
three-dimensional)  vortex  shedding  disappears.  An  additional  factor  tending  to 
complicate  the  experimental  resolution  of  these  issues  is  that  in  either  a  fixed  reference 
frame  or  one  translating  with  the  cylinder,  the  generation  and  shedding  of  vortices  is 
easily  masked  (Perry,  Chong  &  Lim  1982)  at  large  values  of  a  by  the  high  velocities 
induced  in  the  near  wake  by  the  rapidly  rotating  cylinder.  Although  experimental 
techniques  for  measuring  vorticity  are  under  development  for  two-  and  three- 
dimensional  flows  (cf.  Klewicki  &  Falco  1991 ),  the  vorticity  distribution  frequently  can 
be  determined  conveniently  by  direct  computation  of  the  vorticity,  which  ;s  a  frame- 
invariant  quantity. 

To  date,  however,  most  of  the  theoretical  studies  have  shed  no  light  on  the  question 
of  whether  cylinder  rotation  can  suppress  vortex  shedding.  The  analytical  investigations 
of  flow  past  a  rotating  and  translating  circular  cylinder  (Krahn  1955;  Glauert  1 957  a,  b ; 
Moore  1957;  and  Wood  1957)  are  based  on  steady  boundary-layer  theory,  and  are 
hence  inapplicable  to  investigation  of  the  unsteady  separated  flow  associated  with 
vortex  shedding.  The  computational  investigations  of  flow  generated  by  a  rotating  and 
translating  cylinder  reported  by  Ta  Phuoc  Loc  (1975),  Lyulka  (1977),  Townsend 
(1980),  Ingham  (1983),  Ingham  &  Tang  (1990),  and  Tang  &  Ingham  (1991)  concern 
only  the  steady  flow  with  Re  30.  Although  Shkadova  (1982)  discussed  a 
computational  algorithm  for  the  unsteady  flow,  she  presented  only  a  single  set  of 
streamlines  for  each  of  a  few  combinations  of  Re  and  a  (Re  =  20,  40,  and  80  for 
a  =  0.2,  and  Re  =  40  for  a  =  3).  She  did  not  discuss  unsteady  effects  for  the  case  of  a 
rotating  and  translating  cylinder,  and  it  is  not  clear  whether  the  streamlines  presented 
for  each  combination  of  Re  and  a  pertained  to  a  computed  steady  flow,  or  to 
instantaneous  streamlines  (at  unspecified  times)  in  an  unsteady  flow.  In  the  earlier 
work  of  Simuni  (1967)  concerning  the  flow  generated  by  a  cylinder  accelerated 
smoothly  (rather  than  impulsively)  into  rotatory  and  rectilinear  motion,  the  time- 
dependence  of  the  body  motion  was  not  clearly  specified,  nor  was  any  information 
provided  about  the  time-dependence  of  the  computed  solution. 

To  the  best  of  our  knowledge,  other  than  the  experimental  work  of  Coutanceau  & 
Menard  (1985),  Badr  et  al.  (1990),  and  Kimura  et  al.  (1992),  the  only  investigations  of 
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vortex  shedding  by  a  steadily  routing  and  translating  circular  cylinder  in  the  laminar 
regime  are  the  computation:.',  studies  of  Badr  &  Dennis  (1985),  Badr.  Dennis  &  Young 
(1989),  Badr  et  al.  (1990).  Kimura  et  al.  (1992)  and  Chang  &  Chem  (1992).  Although 
Badr  and  coworkers  realized  the  importance  of  extending  the  computations  to  larger 
a,  their  work  for  Re  -  200  was  limited  to  a  ^  1,  and  so  shed  no  light  on  the  conclusion 
of  Coutanceau  &  Menard  (1985)  cited  above.  The  computational  work  of  Kimura 
et  al.,  covering  the  Re  range  from  approximately  400  to  10\  uses  a  discrete  vortex 
method  with  the  cylinder  surface  divided  into  only  14  segments.  The  smallest  Re  at 
which  these  authors  present  results  is  at  a  poorly  defined  value  near  400,  for  which  they 
indicate  only  that  “meandering’  of  the  wake  occurs  for  a  $  1.8,  and  that  for  a  2.6 
‘meandering  disappears,  but  this  value  is  not  so  definitive’.  On  the  other  hand, 
although  the  extensive  and  well-resolved  two-dimensional  computations  of  Chang  & 
Chem  support  the  authors'  detailed  description  of  different  two-dimensional  flow 
regimes  in  the  range  10s  <  Re  <  10*.  a  <  2,  they  concern  a  range  of  Re  and  a  in  which, 
on  the  basis  of  what  is  known  of  the  non-rotating  case  (Williamson  1988).  the  two- 
dimensionality  of  the  flow  is  in  doubt. 

In  the  present  work,  for  Re  =  200,  we  present  computations  extending  the  range  of 
cl  to  the  highest  value  (3.25)  studied  experimentally  by  Coutanceau  &  Menard  (1985). 
After  comparing  our  results  at  lower  values  of  a  to  the  previous  work  of  Coutanceau 
&  Menard  (1985)  and  Badr  &  Dennis  (1985),  we  present  results  for  the  two  largest 
values  of  a  (2.07  and  3.25)  considered  by  the  former  authors.  Ours  is  the  first 
computational  investigation  at  Re  —  200  for  a  >  1,  and  is  significant  in  the  light  of  the 
earlier  conclusions  that  vortex  shedding  is  suppressed  for  a  >  1  or  a  >  2.5.  Our 
computed  results  are  in  excellent  agreement  with  the  experiments  of  Coutanceau  & 
Menard  (1985)  for  Re  =  200  and  a.  —  2.07.  We  then  present  strong  evidence  in  support 
of  the  hypothesis  that  rotation  does  not  suppress  vortex  shedding  for  Re  -  200  and 
a  =  3.25.  This  evidence,  consisting  of  streamlines  viewed  from  a  reference  frame 
moving  with  a  vortex,  and  of  contours  of  constant  vorticity,  is  of  a  type  not 
easily  obtainable  in  the  laboratory,  and  provides  an  important  demonstration 
of  the  capability  of  computational  methods  to  resolve  questions  arising  from 
experiment. 

In  §2,  we  present  the  governing  equations,  along  with  a  transformed  version 
appropriate  for  computations  with  a  body-fitted  time-dependent  grid  used  for 
a  =  3.25.  In  §3,  we  briefly  describe  the  numerical  methods  employed,  including  a  new 
implementation  of  the  Biot-Savart  law  used  to  satisfy  boundary  conditions  on  the 
vorticity  at  the  cylinder  and  on  the  velocity  in  the  far  field.  Section  4  discusses  a 
technique,  more  general  than  that  employed  in  earlier  studies,  for  determining  the 
initial  flow  (at  t  =  0+).  The  main  results,  inluding  comparison  to  the  experiments  of 
Coutanceau  &  Menard  (1985)  and  discussion  of  features  of  the  flow  not  elucidated  by 
their  work  or  earlier  computations,  are  presented  in  §5,  followed  by  a  more  general 
discussion  in  §6. 


2.  Governing  equations 

A  non-rotating  reference  frame  translating  with  the  cylinder  is  used  in  this  study.  In 
this  frame,  the  fluid  at  infinity  has  a  uniform  velocity  of  magnitude  U  in  the  .v-direction, 
and  the  cylinder  rotates  in  the  counterclockwise  direction  with  angular  velocity  C2ez.  as 
shown  in  figure  1 . 

We  use  a  velocity/vorticity  formulation,  consisting  of  the  vorticity  transport 
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Figure  1.  Definition  sketch. 


equation  and  a  vector  Poisson  equation  for  velocity.  In  two  dimensions,  the 
dimensionless  equations  are  (Wu  1975;  Fasel  1976): 

=  (1) 

at  Re 

and  V!F=-CxK),  (2) 

where  we  use  the  cylinder  radius  a  as  the  lengthscale.  and  a/U  as  the  timescale.  The 
velocity  is  normalized  by  U.  Equation  (2)  is  derived  from  the  continuity  equation 

VV=0,  (3) 

the  definition  of  vorticity  for  a  two-dimensional  flow 

(oez  —  V  x  V,  (4) 

and  the  vector  identity 

VxVx  V=  V(V-F)-V2K,  (5) 

where  V  —  uex  +  vey  is  the  velocity  vector. 

The  dimensionless  boundary  conditions  are 

V  =  ex  at  infinity  (6) 

and  V  =  a(  -  ex  sin  6  +  ey  cos  6)  on  the  cylinder  surface.  (7) 

To  allow  for  computation  of  the  flow  on  a  time-dependent  grid,  we  write  (1)  and  the 
components  of  (2)  in  general  body-fitted  (£,  ^-coordinates  as 

=  j  [*,(w£>’,  -  w,  >’f)  -  >v(wr  -X,  -  w,  x,)  -  y,(u«)5 + y£uio\  -  x,(M£ + x^rw),] 

+  (Sa>u  -  2p(i)^  +  y«„)  +  £?%),  (8) 

Sua  -  20U',  +  yum  +  J\Pu,  +  Qun)  =  J((o<  xf  - -  w,  x{),  (9) 

8V'i-20V',+yv„+J\Pvg+Qvlt)  =  /(w£y,-w,y{)  (10) 


and 
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(Reddy  &  Thompson  1977)  where 

P  =  xix'+y\y?  Y  =  xl+>i\  (11) 

P  =  izx+ivv,  Q  =  Vxx  +  Vyv  J 

and  J  =  xiy,-xlly(  (12) 

is  the  Jacobian  of  the  mapping  between  the  (*,>>)-  and  (£,  »/)-coordinate  systems.  Here, 
subscripts  denote  partial  differentiation.  In  (8)  and  in  the  computer  code  developed,  we 
have  allowed  the  grid  in  the  physical  (x,y)- space  to  be  time-dependent.  This  introduces 
additional  terms  associated  with  xt  and  yt  into  the  governing  equations  in  the 
generalized  coordinate  system.  In  this  work,  the  body-fitted  grid  is  simply  one  of 
cylindrical  polar  coordinates  and  is  time-independent,  except  for  a  =  3.25  where  the 
grid  is  made  time-dependent  for  24  ^  <  54.  The  grid  is  uniformly  spaced  in  the 
circumferential  direction  and  is  stretched  in  the  radial  direction,  as  described  below. 


3.  Numerical  methods 

In  this  and  other  numerical  simulations,  it  is  necessary  to  confine  the  computation 
to  a  finite  spatial  domain.  As  a  result,  (6)  cannot  be  applied  directly  at  the  outer 
perimeter  of  the  computational  domain.  Various  far-field  boundary  conditions, 
including  those  derived  from  potential  flow  and  Oseen  expansions,  have  been  adopted 
in  the  past.  The  conditions  imposed  at  the  outer  perimeter  of  the  computational 
domain  have  been  found  to  strongly  influence  the  accuracy  of  steady  flow  computations 
in  this  unbounded  geometry  (Fomberg  1980;  Ingham  1983;  Ingham  &  Tang  1990).  A 
second  difficulty,  common  to  most  simulations  based  on  primitive  variable  (pres¬ 
sure/velocity)  or  velocity/vorticity  formulations,  is  that  conditions  on  either  the 
pressure  or  vorticity  are  required  at  solid  boundaries.  In  this  work,  both  of  these 
difficulties  are  resolved  by  use  of  a  new  implementation  of  the  Biot-Savart  law  briefly 
described  below.  For  further  details,  the  reader  is  referred  to  the  work  of  Wu  and 
coworkers  (Wu  &  Thompson  1973;  Wu  1976;  Wang  &  Wu  1986),  and  Chen  (1989). 

The  definition  (4)  allows  determination  of  the  vorticity  field  from  a  known  velocity 
field.  Conversely,  one  can  determine  the  velocity  field  from  a  known  vorticity  field  via 
the  generalized  Biot-Savart  law,  which  in  two  dimensions  can  be  written  as 


V(r0,t) 


-us. 


w(r,  t)ez  x  (r— r0) 


r-rn 


d  A 


-HI, 


2Q(t)ez  x  (r— r0) 


r-r„ 


dA  +  Vx  (13) 


(Payne  1958;  Wu  1976),  where  the  subscript  0  denotes  the  field  point  where  the  velocity 
is  evaluated,  and  Vx  is  the  uniform  flow  at  infinity.  The  first  double  integral  in  (13)  is 
evaluated  numerically  over  the  fluid  domain  D,  while  the  second  is  evaluated 
analytically  over  the  solid  body  B.  Here,  wez  is  the  vorticity  at  a  point  within  the  fluid, 
and  Qez  is  the  angular  velocity  of  a  point  within  B. 

Equation  (13)  serves  two  purposes  in  this  study.  First,  if  the  vorticity  field  w(r,  t)  is 
known  and  the  domain  D  is  large  enough  to  contain  all  of  the  vorticity  generated  at 
the  solid  boundary  prior  to  time  r,  then  the  velocity  on  the  perimeter  of  D  can  be 
evaluated  directly  by  numerical  integration  of  (13).  Second,  by  linking  the  velocity  and 
vorticity  fields,  (13)  provides  a  basis  for  determining  the  vorticity  on  the  solid 
boundary.  Applying  (13)  to  points  rb  on  the  solid  boundary,  one  obtains 


nr, 


■"-HI. 


w(r.  t)ezx(r-rb) 

lr— rol2 


d  A 


-ill. 


2  a(r)ezx(r-rb) 
lr  —  re>l2 


d A+Vx.  (14) 
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f  Vx  and  the  body  motion  V(r„,t)  and  Q(t)  are  given,  and  the  vorticity  is  known 
:verywhere  except  on  the  solid  boundary,  then  the  only  unknown  in  (14)  is  the 
/Orticity  on  the  solid  boundary.  Therefore,  one  can  solve  (14)  as  a  vector  Fredholm 
ntegral  equation  of  the  first  kind  to  obtain  the  vorticity  values  at  grid  points  on  the 
.olid  boundary. 

In  this  work,  the  numerical  integration  of  (13)  and  (14)  is  performed  over  each 
quadrilateral  element  in  D  using  isoparametric  representations  commonly  used  in 
finite-element  computations.  All  variables  are  located  at  the  intersections  of  grid  lines, 
namely  comers  of  quadrilateral  elements.  The  vorticity  distribution  over  each 
quadrilateral  element  is  approximated  by  bilinear  shape  functions.  Integration  is  then 
performed  over  the  [—  1, 1]  x  [—  1, 1]  square  in  the  isoparametric  plane.  When  the  field 
point  is  far  from  the  quadrilateral  element,  more  efficient  asymptotic  formulae  (Weston 
&  Liu  1982;  Ting  1983)  are  employed.  Further  details  are  provided  in  Chen  (1989). 

In  deriving  the  discretized  forms  of  (8H10),  second-order  central  differences  are 
used  for  all  derivatives  in  the  radial  direction  if,  while  a  pseudospectral  transform 
method  (Orszag  1980;  Zang,  Wong  &  Hussaini  1982)  is  used  to  evaluate  all  derivatives 
in  the  circumferential  direction  £.  The  cross-derivative  terms  are  approximated  by 
central  differencing  in  j?  followed  by  pseudospectral  transformtion  in  £.  We  use  a  fully 
explicit  method  to  advance  the  vorticity  transport  equation  (8)  in  time  to  obtain  the 
vorticity  values  at  the  interior  grid  points.  The  vorticity  on  the  outer  perimeter  of  the 
computational  domain  is  obtained  by  extrapolation.  At  this  stage,  the  vorticity  on  the 
solid  boundary  lags  by  one  time  step.  If  this  vorticity  field  were  used  to  evaluate  the 
right-hand  side  of  (14),  the  result  would  not  satisfy  the  no-slip  and  no-penetration 
conditions.  The  continuous  generation  of  vorticity  at  the  cylinder  is  properly  simulated 
in  our  computations  by  adding  or  subtracting  vorticity  at  the  boundary  at  each  time 
level  in  order  to  satisfy  (14)  identically. 

It  should  be  noted  that  the  solution  of  (14)  is  not  unique  (Wu  1976;  Taslim,  Kinney 
&  Paolino  1984).  To  render  the  solution  unique,  Wu  (1976)  developed  and  imposed  the 
principle  of  vorticity  conservation,  which  states  that  the  total  vorticity  in  the  combined 
fluid  and  solid  regions  must  be  zero  at  all  times.  A  more  general  and  robust  procedure, 
applicable  to  flows  over  single  and  multiple  solid  bodies,  has  been  developed  by  Chen 
(1989),  and  is  used  here. 

The  computational  loop  to  advance  the  solution  from  one  time  level  to  the  next 
consists  of  the  following  steps; 

(i)  The  discretized  vorticity  transport  equation  (8)  is  advanced  explicitly  to  obtain 
new  vorticity  values  at  all  interior  grid  points,  using  a  second-order  rational 
Runge-Kutta  method  (Wambecq  1978).  In  contrast  to  the  three-step  Adams- 
Bashforth  method  used  by  previous  authors,  this  allows  a  much  larger  time-step  size 
due  to  a  less  severe  stability  constraint. 

(ii)  Using  known  vorticity  values  at  the  interior  grid  points,  the  kinematic  constraint 
(14)  is  used  to  update  the  vorticity  values  on  the  solid  boundary. 

(iii)  Using  the  updated  vorticity  field,  the  velocity  at  points  on  the  outer  perimeter 
of  the  computational  domain  is  evafuated  from  (13).  Once  velocity  boundary  values  are 
known,  Poisson  equations  (9)  and  (10)  are  solved  for  the  new  velocity  field.  The 
discretizations  of  (9)  and  (10)  are  11 -banded  and  block-diagonal  in  form,  and  are 
solved  by  a  preconditioned  biconjugate  gradient  algorithm  (Chen,  Koniges  &  Anderson 
1989). 

The  method  outlined  above  is  particularly  well-suited  for  the  initial  development  of 
the  flow  generated  by  impulsively  started  bodies.  This  is  so  because  the  vorticity  is 
initially  concentrated  near  the  solid  body,  thus  allowing  the  numerical  simulation  to  be 
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confined  to  a  domain  containing  nearly  all  the  vorticity.  Since  the  velocity  at  the  outer 
perimeter  of  the  computational  domain  is  calculated  via  (13),  imposition  of 
computational  far-field  velocity  boundary  conditions  is  avoided.  The  use  of  (13)  in 
calculating  the  far-field  velocity  would  be  numerically  exact  at  time  t  if  the 
computational  domain  contained  all  the  vorticity  generated  at  the  cylinder  surface 
prior  to  t.  We  note  that,  in  the  streamfunction/vorticity  formulation,  the  stream- 
function  values  on  the  outer  perimeter  of  the  computational  domain  can  be 
obtained  similarly  by  using  an  integral  constraint  equation  (Wu  &  Sampath  1976).  We 
also  note  that  the  velocity  field  can  be  obtained  by  applying  (13)  at  every  grid  point. 
This  results  in  a  point-by-point  scheme  which,  unlike  schemes  using  the  Poisson 
equations  (9)  and  (10),  does  not  require  solution  of  large  linear  equation  systems.  This 
approach  was  adopted  in  Wu’s  earlier  work,  and  might  be  attractive  on  massively 
parallel  computers. 

The  size  of  the  computational  domain  is  chosen  according  to  the  time  span 
investigated.  Here,  we  set  the  outer  boundary  of  the  computational  domain  to  be  a 
circle  of  radius  24  for  t  ^  24.  For  a  =  3.25.  the  grid  is  made  time-dependent  for 
24  ^  t  ^  54.  We  use  128  uniformly  spaced  and  120  stretched  grid  lines  in  the  ft-  and  re¬ 
directions,  respectively.  The  stretching  function  of  Vinokur  (1983)  is  used  to  distribute 
the  circular  grid  lines  on  1  ^  r  <  24.  This  stretching  allows  grid  points  to  be  clustered 
near  one  or  both  ends  of  the  domain,  or  anywhere  between,  by  adjusting  two 
parameters  s0  and  sv  Here,  s0  and  s1  are  the  ratios  of  the  spacing  if  N  points  were 
distributed  uniformly,  to  the  actual  spacings  at  the  inner  and  outer  boundaries, 
respectively.  For  a  =  0.5  and  1.0,  we  set  s0  =  5.0  and  r,  =  0.25.  The  grid  spacings 
adjacent  to  the  cylinder  and  at  the  outer  perimeter  are  4%  and  75%  of  the  cylinder 
radius,  respectively.  For  a  —  2.07  and  3.25,  we  set  sa  =  10.0  and  s,  =  0.25  to  further 
cluster  circular  grid  lines  near  the  cylinder,  with  grid  spacings  adjacent  to  the  cylinder 
and  at  the  outer  perimeter  being  2%  and  76%  of  the  cylinder  radius,  respectively.  At 
the  end  of  the  simulation,  the  vorticity  magnitude  on  the  outer  perimeter  is  less  than 
10~5  for  all  cases  except  as  noted,  indicating  that  only  a  negligible  fraction  of  the 
vorticity  has  escaped  the  domain. 

With  the  grid  chosen  for  a  $  2.07,  our  method  requires  approximately  5  CPU 
seconds  per  time  step  on  a  CRAY  2. 


4.  Determination  of  the  initial  flow 

In  most  numerical  simulations  of  flow  over  impulsively  started  bodies,  the  initial 
flow  field  is  taken  to  be  the  potential  flow,  since  the  vorticity  at  t  =  0+  is  concentrated 
on  the  body  surface  in  the  form  of  a  vortex  sheet.  Perturbation  solutions  in  which  t  is 
the  small  parameter  have  also  been  used  as  initial  conditions  (Collins  &  Dennis  1973 
for  a  =  0;  Badr  &  Dennis  1985  for  a  +  0).  Both  approaches  require  special  techniques 
in  order  to  obtain  the  initial  flow  field.  In  the  present  work,  determination  of  the  initial 
flow  field  requires  no  special  treatment.  The  same  procedure  used  for  determining  the 
boundary  vorticity  distribution  satisfying  the  no-slip  and  no-penetration  conditions  is 
applied.  More  specifically,  (14)  is  solved  for  the  unknown  boundary  vorticity  at  t  =  0+, 
with  the  vorticity  taken  to  be  zero  at  every  grid  point  away  from  the  cylinder  surface. 
Once  the  boundary  vorticity  values  are  obtained,  the  initial  velocity  field  is  determined 
by  solving  (9)  and  (10),  with  the  velocity  on  the  outer  perimeter  of  the  computational 
domain  determined  by  application  of  (13)  to  points  on  the  outer  perimeter.  This 
versatile  procedure  enables  the  numerical  code  to  handle  bodies  of  arbitrary  shape 
undergoing  arbitrary  rotational  and  translational  motion. 
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Errors  in  approximating  the  vorticity  layer  of  infinitesimal  thickness  at  /  =  0*  are 
inherent  to  any  computational  scheme.  However,  as  discussed  by  Lugt  &  Haussling 
(1974),  even  in  the  worst  case  of  a  body  set  impulsively  into  motion,  the  duration  of 
these  errors  is  confined  to  a  very  limited  time  close  to  /  =  0.  For  the  present  algorithm, 
this  will  be  confirmed  in  §5  by  comparison  at  small  times  of  our  numerical  results  to 
the  perturbation  solution  of  Badr  &  Dennis  (1985). 


5.  Results 

In  this  section,  numerical  results  for  Re  =  200  with  a  =  0.5,  1.0.  2.07,  and  3.25  are 
presented  and  discussed.  The  parameter  values  chosen  allow  comparison  to  the 
experimental  results  of  Coutanceau  &  Menard  (1985)  and  permit  a  critical  evaluation 
of  their  conclusions  regarding  the  suppression  of  vortex  shedding  by  rotation. 

For1  a  -  0.5  and  1.0,  excellent  agreement  of  our  computed  streamlines  with  previous 
experimental  (Coutanceau  &  Menard  1985)  and  numerical  (Badr  &  Dennis  1985) 
results  is  obtained.  For  a  =  2.07,  for  which  no  numerical  results  have  been  reported 
previously,  we  obtain  excellent  agreement  with  experiment.  For  a  =  3.25,  the  relatively 
small  disagreement  between  experiment  and  our  computations  is  probably  due  to  the 
effects  of  three-dimensionality  and  sidewall  boundary  layers  in  the  former  (see  §5.4). 
For  a  =  0.5,  1.0,  and  2.07,  we  continue  the  simulations  to  larger  dimensionless  times 
(t  —  24)  than  could  be  studied  experimentally  by  Coutanceau  &  Menard  (1985),  so  that 
the  nature  of  the  wake  development  can  be  better  discerned.  For  a  =  3.25,  we  extend 
our  calculation  to  t  =  54  to  include  shedding  of  the  second  and  third  vortices. 

For  the  same  values  of  a,  we  also  present  computations  of  the  vorticity  contours, 
which  we  find  useful  for  studying  vortex  shedding  without  the  effects  of  ‘masking’ 
associated  with  frame-dependent  streamlines.  We  also  present  trajectories  of  the  shed 
vortices,  computed  from  vorticity  contours  and  from  streamlines. 

We  note  that  the  timescale  adopted  here  is  the  same  as  that  used  by  Badr  &  Dennis 
(1985).  Conversion  of  the  dimensionless  times  of  Coutanceau  &  Menard  (1985)  to 
those  herein  requires  multiplication  of  the  former  by  a  factor  of  two. 

As  discussed  in  §4,  the  solution  procedure  presented  in  §3  is  applied  at  t  =  (T  to 
obtain  the  initial  flow  field.  Errors  are  present  due  to  the  inability  of  any  numerical 
scheme  to  resolve  the  infinitesimal  vorticity  layer  at  /  =  0+.  To  confine  these  errors  to 
small  times  near  t  =  0,  small  initial  time  steps  are  used  for  all  cases.  For  the  first  20  time 
steps.  At  —  10-4  is  used.  This  is  followed  by  28  steps  with  At  =  10-3,  which  brings  the 
time  level  in  0.03.  A  constant  At  (10_i  for  a  =  0.5  and  1.0;  2.5  x  10-3  for  a  =  2.07  and 
3.25)  is  taken  for  the  rest  of  each  simulation.  The  time-step  size  is  not  dictated  by 
the  numerical  stability  constraint,  but  rather  is  chosen  on  the  basis  of  accuracy 
considerations.  To  demonstrate  the  accuracy  of  the  initial  flow  field,  we  show  in  figure 
2  (a-d)  the  variation  of  vorticity  on  the  cylinder  for  four  values  of  a  at  small  times.  The 
numerical  results  are  compared  to  the  asymptotic  formula 

w(  1,0, /)  *  ^|a^A-0  +  |A  +  Ajsin0+^2.7844A— ^Aja/cos0 

/sin  2^  (15) 


given  by  Badr  &  Dennis  (1985),  where 

A  =  (8  t/ Re}. 


■  6-55™-5(,4) 


(16) 
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Figure  2  (a,  b).  For  caption  sec  facing  page. 


We  see  that  agreement  with  the  asymptotic  results  improves  as  t  increases.  Better 
agreement  is  achieved  for  a  =  2.07  and  3.25  than  for  smaller  a  since  a  smaller  time-step 
size  is  used.  These  results  demonstrate  that  errors  are  indeed  confined  to  a  very  limited 
time  close  to  t  =  0,  as  reported  by  Lugt  &  Haussling  (1974). 

5.1.  Results  for  a  =  0.5 

In  this  subsection,  we  extend  the  a  =  0.5  computations  of  Badr  &  Dennis  (1985)  to 
larger  dimensionless  times  than  considered  by  them  or  in  the  experiments  of 
Coutanceau  &  Menard.  The  kinematics  and  dynamics  of  vortex  shedding  are  discussed 
using  instantaneous  streamlines  in  two  different  reference  frames,  as  well  as  vorticn., 
distributions.  The  streamfunction  is  computed  from  the  velocity  field  by  a  least-squares 
method  described  by  Chen  (1989). 

Computations  performed  for  0  ^  t  ^  24  show  that  the  results  are  in  excellent 
agreement  with  the  experimental  work  of  Coustanceau  &  Menard  (1985)  for 
1  <  /  <  13  and  the  computations  of  Badr  &  Dennis  (1985)  for  1  ^  t  ^  12. 
Figure  3  (a-j)  shows  a  sequence  of  instantaneous  computed  streamlines  for  8  <  /  <  24, 
viewed  from  a  non-rotating  frame  translating  with  the  circular  cylinder.  We  adopt  the 
notation  used  by  Coutanceau  &  Menard  (1985)  in  their  discussion  of  the  flow 
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Figure  2.  Evolution  of  the  vorticity  distribution  on  the  surface  of  the  cylinder  at  early  times. 

Re  =  200.  Symbols:  asymptotic  solution;  □,  /  =  0.05;  A,  t  =  0.1 ;  O,  t  =  0.2.  - ,  Numerical 

solution,  (a)  a.  —  0.5,  (6)  a  =  1.0,  (c)  a  =  2.07,  ( d )  a  =  3.25. 


development  for  t  ^  13,  to  which  the  reader  is  referred.  We  observe  that  the 
coalescence  of  two  ‘intermediate’  eddies  E,  and  E3  to  form  eddy  E3  at  t  %  12  (figure 
3  c,  d)  indeed  repeats  (at  t  &  22;  figure  3  h,  i )  as  predicted  by  Coutanceau  &  Menard 
(1985).  (The  subscript  here  denotes  the  order  of  appearance  of  eddies  after  the 
impulsive  start.)  The  transposition  of  saddle  points  S2  and  S3  associated  with  eddies  E, 
and  E'3,  respectively,  discussed  by  Coutanceau  &  Menard  (1985)  and  sketched  by  Badr 
et  al.  (1986),  is  clearly  shown  in  figure  3  (a,  b).  Also,  a  common  boundary  for  E3  and 
E3,  which  was  difficult  to  observe  experimentally  (Coutanceau  &  Menard  1985)  due  to 
limitations  of  the  flow  visualization  technique,  does  indeed  exist,  as  shown  in  figure 
3(a).  The  common  boundary  soon  becomes  an  ‘alleyway’  for  fluid  to  pass  through,  as 
shown  in  figure  3(6).  As  noted  by  Eaton  (1987),  existence  of  such  an  alleyway  in  an 
unsteady  flow  does  not  imply  that  fluid  is  carried  from  one  side  of  the  wake  to  the 
other. 

As  discussed  by  Perry  et  al.  (1982),  the  streamlines  are  not  invariant  with  respect  to 
a  change  in  reference  frame,  and  the  vortex  street  can  appear  very  differently  in 
different  frames.  To  best  observe  the  development  of  a  vortex,  the  observer  should 


move  with  its  centre  (Lugt  1979).  Otherwise,  the  vortex  can  be  masked  by  the  motion 
of  the  observer  relative  to  the  vortex  (Williamson  1985;  Coutanceau  &  Menard  1985). 
This  masking  phenomenon  was  described  by  Coutanceau  &  Menard  as  the  opening  up 
of  vortices  and  disappearance  of  closure  points  into  a  wave-like  pattern,  as  sketched 
in  Lugt  (1979).  Since  an  attached  vortex  translates  with  the  cylinder,  it  can  be  clearly 
observed  in  a  frame  translating  with  the  cylinder.  However,  after  the  vortex  is  shed,  its 
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(0  .v.:  .1: -  :  ;;;  U) 


Figure  3.  Instantaneous  streamlines  for  Re  =  200,  at  =  0.5  at  various  times,  viewed  from  a  non¬ 
rotating  frame  translating  with  the  cylinder.  Streamlines  with  non-negative  (including  zero)  and 
negative  streamfunction  values  (v>)  are  shown  by  solid  and  dashed  lines,  respectively.  At  each  time, 
the  cylinder  is  a  streamline  with  ijr  *  0.  The  values  of  ^  plotted  are  0,  —0.01,  ±0.02,  —0.03.  ±0.04. 
±0.06.  ±0.08.  ±0.10,  ±0.12,  ±0.15,  ±0.20.  ±0.25,  ±0.30,  ±0.35,  ±0.40,  ±0.45,  ±0.50.  ±0.60. 
±0.70.  ±0.80.  ±  1.00,  with  an  increment  of  ±0.2  thereafter,  (a)  r  =  8.0,  (b)  9.0,  (c)  11.0.  (d)  12.0, 
(e)  14.0,  (/)  16.0,  (g)  18.0,  (h)  20.0,  (0  22.0,  (y)  24.0. 

core  is,  especially  in  the  far  wake,  essentially  stationary  with  respect  to  the  free  stream. 
Therefore,  it  is  generally  easier  to  observe  shed  vortices  in  a  frame  fixed  with  the 
undisturbed  fluid.  The  instantaneous  streamlines  observed  in  such  a  frame  are  shown 
at  selected  times  in  figure  4  (a~c).  As  expected,  the  shed  vortices  are  clearly 
distinguishable.  We  note  that  two  additional  vortices  are  shed  over  an  interval  of  about 
10  dimensionless  time  units,  and  that  the  flow  near  the  cylinder  is  very  similar  in  figures 
4(a)  and  4(c).  We  further  note  that  in  a  moving  frame,  the  cylinder  itself  is  not  a 
streamline,  and  the  attached  vortices  in  figure  4(a-c)  are  now  masked  by  the  velocity 
field  in  the  near  wake  of  the  cylinder.  In  a  frame  translating  with  the  cylinder,  however, 
the  shed  vortices  are  hidden  in  the  oscillating  wake,  as  shown  in  figure  5(a-c)  for  the 
values  of  t  shown  in  figure  4(o-c). 

At  the  same  dimensionless  times,  figure  6(a-c)  shows  the  corresponding  vorticity 
contours,  which  are  invariant  with  respect  to  translation  of  the  observer.  Because 
rotation  divides  the  surface  of  the  cylinder  into  ‘downstream-moving’  (n  <6  <  2n) 
and  ‘upstream-moving’  (0  <  6  <  n)  parts,  a  basic  symmetry  of  the  vorticity  field  for 
the  non-rotating  (a  =  0)  case  (associated  with  the  fact  that  for  a  T-periodic  flow,  the 

relations  u (r>8,t)  =  u(r,  -d,t-T),  v[r,d.  t)  =  ~v{r,  -d.t-r)  (17) 

lead  to  w(r.8,r)  =  —<o(r,  —8,t—r),  where  0  <  r  <  T  is  a  phase  difference)  is  broken. 
Nonetheless,  the  process  is  topologically  similar  to  the  a  =  0  case,  with  the  saedding 
of  vortices  of  alternating  rotational  sense  being  associated  with  the  thinning  and 
severance  of  elongated  vorticity  contours  emanating  from  opposite  sides  of  the 
cylinder.  As  in  the  experiments  of  Diaz  et  at.  (1983)  at  higher  Re.  cylinder  rotation  has 
the  effect  of  altering  the  initial  trajectories  of  the  shed  vortices,  although  these  are 
expected  to  become  parallel  to  the  direction  of  cylinder  translation  as  the  vortices  move 
farther  away  from  the  rotating  cylinder  and  are  adverted  downstream.  As  in  the 
a  =  0  case,  vortices  of  opposite  sense  lie  on  opposite  sides  of  a  ‘  street  ’,  although  the 
rotation  has  clearly  displaced  the  midline  of  the  street  upward. 

Figure  7  (a,  b)  shows  that  the  computation  of  vortex  core  and  saddle  trajectories 
using  streamfunction  values  (in  excellent  agreement  with  the  trajectories  computed  by 
Coutanceau  &  Menard  1985  from  experimental  streamline  data  and  shown  in  their 
figure  4a)  can  differ  significantly  from  those  computed  using  the  vorticity  distribution. 


Figure  4.  Instantaneous  streamlines  for  Re  =  200,  a  =  0.5,  viewed  from  a  frame  fixed  with  the 
undisturbed  fluid.  Dashed  (solid)  lines  represent  constant  non-negative  (negative)  streamfunction 
values  with  increments  of  =  ±0.1,  including  ft  =  0.  (a)  t  =  12.0,  (6)  17.0,  (c)  22.0. 

For  example,  at  /  =  7,  the  stream  wise  location  of  the  core  of  the  first  vortex  is  at  about 
x/a  =  3.3  as  determined  from  the  streamfunction  (in  a  frame  translating  with  the 
cylinder),  and  a  bit  less  than  x/a  =  2.5  as  determined  from  the  vorticity.  This  clearly 
illustrates  the  effect  of  streamline  masking  on  vortices  moving  with  velocities 
significantly  different  than  that  of  the  frame  to  which  the  motion  is  referred. 


Figure  5.  Instantaneous  streamlines  for  Re  =  200,  a  =  0.5,  viewed  from  the  reference  frame 
described  in  figure  3.  The  plotting  convention  and  contour  levels  are  as  in  figure  3.  (a)  t  =  12.0,  (b) 
17.0,  (c)  22.0. 
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The  trajectories  of  figure  7  (b)  also  show  that  the  vortex  cores  execute  motions  much 
more  complicated  than  would  be  inferred  from  the  trajectories  computed  from 
streamfunction  values  (figure  la  of  the  present  work  and  figure  4 a  of  Badr  &  Dennis 
1985).  Specifically,  for  the  first  vortex  core,  figure  1(b)  shows  that  although  the  v- 
component  of  the  core  velocity  increases  to  nearly  the  free-stream  value  as  the  vortex 
moves  farther  behind  the  cylinder,  the  i-component  oscillates  (being  sometimes 
negative)  about  a  decidedly  non-zero  mean,  to  a  distance  at  least  15  cylinder  radii 
downstream. 

As  a  further  check  on  the  accuracy  of  our  results,  we  show  in  figure  8  the  temporal 
evolution  of  profiles  of  the  .v-component  of  the  velocity  along  the  v-axis  below  the 
cylinder  for  t  ^  8.  Very  good  quantitative  agreement  with  experimental  data  taken 
from  Coutanceau  &  Menard  (1985)  is  obtained  near  the  cylinder.  Farther  away,  there 
is  slightly  more  scatter  in  the  experimental  data.  Our  results  along  the  positive  y-axis 
are  graphically  indistinguishable  from  those  of  Badr  &  Dennis  (1985).  obtained  using 
different  numerical  methods. 


5.2.  Results  for  a.  =  1.0 

For  at  =  1.0,  we  have  computed  streamlines  analogous  to  those  of  figure  3  (u-j)  for 
at  =  0.5.  Our  results  are  indistinguishable  from  those  of  Badr  &  Dennis  (1985)  for  the 
range  of  f(l  ^  t  ^  12)  covered  in  their  work.  For  14  ^  t  ^  24,  figure  9 (a-f)  shows 
instantaneous  streamlines  viewed  from  a  non-rotating  frame  translating  with  the 
cylinder,  beginning  with  the  largest  value  of  t  considered  in  the  experimental  work  of 
Coutanceau  &  Menard  (1985).  Unlike  the  at  =  0.5  case  where  the  second  eddy  E2 
appears  at  /  ft  2.0,  here  E2  and  the  third  intermediate  eddy  E'  form  almost 
simultaneously  at  t  ft  6.5,  as  shown  in  the  earlier  experimental  and  computational 
work.  During  the  next  cycle  of  vortex  formation,  however,  E4  appears  before  E':  is 
formed,  as  seen  in  figure  9 (b.  c).  In  general,  the  increase  in  at  tends  to  inhibit  the 
formation  of  the  vortex  at  the  downstream-moving  side  of  the  cylinder,  as  reported  in 
previous  experiments  (Diaz  et  al.  1983;  Coutanceau  &  Menard  1985). 

Figure  10(o-e)  shows  that  the  trajectories  of  the  shed  vortices  for  at  =  1.0  are 
qualitatively  similar  to  those  for  at  =  0.5,  except  that  the  vortices  shed  from  the 
downstream-moving  side  now  lie  above  the  midline  of  symmetry  (8  =  0).  and  due  to 
the  counterclockwise  fluid  motion  generated  near  the  cylinder  by  its  rotation,  will 
remain  above  the  midline  during  their  subsequent  advection  downstream.  Otherwise, 
the  topology  of  the  shedding  process  is  altered  relatively  little,  with  vortices  of 
alternating  sense  being  shed  from  opposite  sides  of  the  cylinder,  and  subsequently- 
being  found  on  opposite  sides  of  an.  albeit  distorted,  ‘street’  as  they  are  advected 
downstream. 

We  have  also  computed  the  temporal  evolution  of  profiles  of  the  v-  and  y- 
components  of  the  velocity  along  the  .v-axis  in  the  wake  of  the  cylinder  for  t  ^  8.  Again, 
good  agreement  with  the  experimental  results  of  Coutanceau  &  Menard  (1985)  is 
obtained. 


5.3.  Results  for  a  =  2.07 

As  at  increases,  the  vorticity  layer  generated  at  the  upstream-moving  side  of  the  cylinder 
intensifies,  resulting  in  even  larger  radial  derivatives.  Consequently,  it  becomes  more 
difficult  to  maintain  accuracy,  as  pointed  out  by  Badr  &  Dennis  (1985).  We  are  able 
to  achieve  accurate  numerical  results  by  using  finer  radial  grid  spacings  near  the 
cylinder,  as  discussed  in  §3.  and  a  smaller  time-step  size,  as  discussed  at  the  beginning 
of  this  section.  Figure  1 1  (a-h)  shows  instantaneous  streamlines  in  the  near  wake  for 
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Figure  6.  Equivorticity  contours  for  Re  —  200.  x  —  0.5.  Dashed  (solid)  lines  represent  constant 
positive  (negative)  vorticity  values  with  a  constant  increment  of  ±0.5.  with  the  magnitude  of  the 
weakest  contour  level  shown  being  0.5.  (a)  t  =  12.0.  ( b )  17.0.  ( c )  22.0. 

3  sg  t  <  24,  viewed  from  a  non-rotating  frame  translating  with  the  cylinder.  As  seen  in 
figure  12  (a,  b),  our  results  are  virtually  identical  to  the  flow  visualizations  of 
Coutanceau  &  Menard  (1985)  for  t  ^  9,  the  longest  dimensionless  time  for  which 
experimental  results  are  available.  The  sequence  of  figures  13  (a-d)  for  t  -  2, 4.  6.  and 
8  shows  that  as  t  increases  the  computed  flow  is  also  in  excellent  agreement  with 
experiment  upstream  of  the  cylinder,  and  becomes  increasingly  less  symmetric  about 
the  .T-axis. 

The  excellent  agreement  between  our  two-dimensional  computations  and  the  flow 
visualizations  of  Coutanceau  &  Menard  (1985)  in  a  single  spanwise  plane  strongly 
suggests  that  for  t  <  9  the  flow  is  quite  two-dimensional,  at  least  near  the  centre  of  the 
span  of  the  cylinder,  and  that  sidewall  boundary-layer  effects  are  unimportant. 

To  better  elucidate  the  vortex  shedding  process,  we  show  in  figure  1 4(a~f)  the 
vorticity  contours  for  t  ^  24.  As  for  smaller  values  of  a,  vortex  shedding  still  occurs  for 
a  =  2.07,  with  vortices  of  alternating  rotational  sense  shed  from  opposite  sides  of  the 
cylinder.  However,  at  this  larger  rotation  rate  the  asymmetry  of  the  process  is  clear. 
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Figure  7.  Trajectories  of  the  cores  C,  and  closure  points  S,  for  Re  *  200.  a  =  0.5  obtained  from 
(«)  instantaneous  streamlines,  and  (b)  equivorticity  contours. 


and  manifests  itself  in  the  considerably  reduced  strength  of  the  vortices  shed  from  the 
downstream-moving  side  of  the  cylinder  (relative  to  those  for  a  =  0.5  shown  in  figure 
6 a-c,  and  relative  to  those  of  opposite  rotational  sense  for  a  =  2.07),  as  well  as  in  the 
fact  that  the  shed  vortices  seem  to  be  forming  a  *  single  file’  line,  rather  than  pairing  off 
(according  to  rotational  sense)  on  opposite  sides  of  a  ‘street’.  The  vortices  shed  from 
the  downstream-moving  side  of  the  cylinder  are  relatively  weak  because  immediately 
after  the  impulsive  start,  this  part  of  the  rotating  boundary  travels  at  about  the  same 
speed  as  the  adjacent  fluid.  The  weak  vorticity  layer  generated  near  6  =  |jc  (shown  in 
figure  2  c)  is  responsible  for  the  weakness  of  the  shed  vortex. 

Figure  15  shows  that  the  trajectory  of  the  vortex  core  C,  determined  from  equi¬ 
vorticity  contours  is  again  very  different  from  the  trajectory  determined  by  Coutanceau 
&  Menard  from  streamline  data  (their  figure  14).  In  particular,  the  velocity  of  the  first 
vortex  core  still  has  a  considerable  ^-component  until  at  least  t  =  24.  This  disagreement 
may  be  due  to  either  the  fact  that  the  earlier  determination  of  the  vortex  core  trajectory 
was  made  from  streamline  data  (discussed  in  §5.1),  or  to  the  ‘confining  wall  effect’  in 
the  experimental  work  (discussed  by  Coutanceau  &  Menard  in  conjunction  with  their 
results  for  a  =  3.25). 
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Figure  8.  Temporal  evolution  of  u- velocity  profiles  (in  the  frame  described  in  figure  3)  along  the 

y-axis  in  the  cylinder  wake  for  Re  *  200,  a  =  0.5.  6  =  fir.  - ,  Numerical  solution.  Symbols. 

experimental  data  of  Coutanceau  &  Menard  (1985);  ▼,  l  =  0.5;  A,  4.0;  if,  5.0;  O,  6.0. 


5.4.  Results  for  a  =  3.25 

For  a  =  3.25.  figure  16 (a-d)  shows  instantaneous  streamlines  in  the  near  wake  for 
t  =  5,  9,  16,  and  24,  viewed  from  a  non-rotating  frame  translating  with  the  cylinder.  In 
the  same  reference  frame,  figure  \l{a-d)  shows  comparisons  of  the  computed 
streamlines  to  unpublished  experimental  results  of  Coutanceau  &  Menard  centred 
farther  upstream.  We  note  that  the  agreement  (compare  also  figure  166  to  figure  1 1  d 
of  Coutanceau  &  Menard  1985),  is  very  good.  However,  at  /  =  9  (the  largest 
dimensionless  time  for  which  Coutanceau  &  Menard  reported  results)  the  com¬ 
putational  results  differ  perceptibly  from  the  flow  visualization  of  Coutanceau  & 
Menard  (1985).  We  believe  that  the  differences  result  from  three-dimensional  and 
sidewall  boundary-layer  effects  in  the  experiments  for  large  values  of  a  at  large  time, 
as  discussed  by  Coutanceau  &  Menard  (1985).  Although  for  a  =  0  the  nominally  two- 
dimensional  flow  at  Re  —  200  is  unstable  with  respect  to  three-dimensional 
disturbances  (Williamson  1988);  no  information  regarding  the  effect  of  rotation  on 
stability  appears  to  be  available. 

Computed  vorticity  contours  are  shown  for  8  sj  t  54  in  figure  18(0-6).  At  /  =  24, 
figure  18(c)  shows  that  the  elongated  vorticity  contour  has  not  yet  been  severed.  To 
investigate  whether  the  vortex  shedding  process  continues  to  larger  t  for  this  value  of 
a,  we  let  the  size  of  the  computational  domain  grow  linearly  in  time  until  t  =  54  in 
order  to  extend  the  computation.  The  same  values  of  s0  and  j,  in  the  stretching  function 
are  used  to  distribute  120  circular  grid  lines  between  r=  1  and  r=  1.5/-  12.  The 
numerical  solution  at  /  =  54  is  not  as  accurate  as  at  smaller  dimensionless  times  since 
the  vorticity  at  the  outer  perimeter  ( r  =  69)  is  not  negligible  (on  the  order  of  10"1  or 
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Figure  9.  Instantaneous  streamlines  for  Re  -  200.  a=  1.0  at  various  times,  viewed  from  the 
reference  frame  described  in  figure  3.  The  plotting  convention  and  streamfunction  values  are  as  in 
figure  3.  (<?)  I  =  14.0.  (b)  16.5.  (c)  17.5.  (d)  18.0.  (e)  21.0,  (/)  24.0. 


smaller)  at  the  end  of  the  simulation.  Nonetheless,  we  believe  that  our  computation  of 
the  vortex  shedding  in  the  near  wake  is  qualitatively  correct,  even  for  t  >  32  (at  which 
time  the  maximum  of  the  absolute  value  of  the  vorticity  on  the  outer  boundary  is  still 
less  than  4  x  10-4). 

Figure  18(e)  shows  that  a  second  vortex  is  shed  into  the  wake  as  for  smaller  a, 
although  shedding  here  occurs  at  a  much  later  time.  However,  even  in  a  reference  frame 
fixed  with  the  undisturbed  fluid,  the  weak  second  vortex  is  masked  by  the  high  velocity 
induced  in  the  near  wake  by  the  rapidly  rotating  cylinder.  Therefore,  for  large  a  it  is 
not  surprising  that  previous  flow  visualization  experiments,  including  those  performed 
in  a  frame  translating  with  the  cylinder,  have  failed  to  reveal  the  shedding  of  a  second 
(weaker)  vortex.  Moreover,  in  the  work  of  Coutanceau  &  Menard  (1985).  experimental 
limitations  prevented  continuation  of  the  flow  to  the  dimensionless  time  at  which  the 
second  vortex  would  have  been  shed.  These  difficulties  have  led  to  the  erroneous 
conclusions  that  for  a  >  2.5  (or  a  >  1).  no  vortices  are  shed  after  the  initial  (strong) 
starting  vortex,  and  that  the  vortex  street  is  completely  destroyed. 
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Figure  10.  Equivorticity  contours  for  Re  —  200,  a  =  1.0.  Dashed  (solid)  lines  represent  constant 
positive  (negative)  vorticity  values  with  a  constant  increment  of  0.5.  with  the  magnitude  of  the 
weakest  contour  level  shown  being  0.5.  (a)  t  =  8.0,  (b)  10.0.  (c)  14.0.  id)  20.0.  (e)  24.0. 


To  render  the  second  vortex  distinguishable  using  streamlines,  we  show  in  figure 
19(a)  the  streamlines  at  t  =  32  in  a  reference  frame  moving  with  the  core  of  the  second 
vortex.  Experimentally,  this  would  be  a  difficult  task  since  a  camera  would  have  to 
move  with  the  vortex  core  velocity  vector  (ueez+veev)  relative  to  the  cylinder,  which 
is  not  known  a  priori.  Numerically,  this  is  achieved  by  simply  subtracting 
streamfunction  values  corresponding  to  the  velocity  of  the  vortex  core  from  those  in 
an  inertial  reference  frame  fixed  with  the  cylinder.  The  third  vortex  is  clearly  visible  in 


Figure  1 1 .  Instantaneous  streamlines  for  Re  =  200,  a  =  2.07  at  various  times,  viewed  from  the 
reference  frame  described  in  figure  3.  In  addition  to  those  contour  levels  shown  in  figure  3. 
t  =  -0.17,  -0.19,  -0.21,  -0,22,  -0.23.  -0.24  are  also  plotted,  (a)  t  =  3.0,  ( b )  7.0,  (c)  13.0,  (d) 
17.0.  (e)  21.0,  (/)  22.0,  (g)  22.5,  (h)  24.0. 

a  frame  translating  with  its  core  (figure  19  b),  while  the  first  vortex  is  still  sufficiently 
strong  to  be  discernible  in  this  reference  frame  (in  which  its  core  velocity  is  non-zero). 

In  contrast  to  the  results  presented  above  for  smaller  a,  for  a  =  3.25  the  third  vortex 
shed  is  of  the  same  rotational  sense  as  the  second.  This  differs  from  the  case  Re  =  103 
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Figure  12.  Comparison  of  computed  (left)  and  experimental  (right)  instantaneous  streamlines  for 
Re  =  200.  3.  =  2.07.  The  camera  in  the  experiment  and  the  reference  frame  in  the  computation 
ranslate  with  the  cylinder,  (a)  i  =  5.0,  (h)  t  =  9.0. 


and  a  =  3  studied  by  Badr  et  al.  (1990),  in  which  the  first  two  vortices  formed  are  of 
the  same  sense,  but  only  one  is  shed.  Moreover,  in  the  two-dimensional  computations 
of  Badr  et  al.  (1990),  only  two  vortices  were  formed,  and  the  computed  flow- 
approached  a  steady  state.  With  respect  to  understanding  the  real  flow,  however,  their 
results  must  be  regarded  with  caution,  as  the  experiments  of  the  same  authors  clearly 
show  that  the  flow  is  three-dimensional  for  Re  =  103  and  a  =  3. 

For  a  =  3.25  the  trajectory  in  figure  20  shows  that,  in  the  mean,  the  first  vortex 
indeed  continues  to  move  upwards,  in  contrast  to  the  result  of  Coutanceau  &  Menard 
(their  figure  14)  in  which  the  first  vortex  apparently  moves  back  to  the  midline 
( y  =  0).  Coutanceau  &  Menard  ascribe  this  artifact  to  the  effect  of  a  confining  (upper) 
wall,  a  complication  not  present  in  our  computations. 

5.5.  Temporal  evolution  of  the  lift  and  drag  coefficients 
Finally,  we  present  the  temporal  evolution  of  the  lift  and  drag  coefficients  defined  by 

CL  =  L/pUia  (18) 

and  CD  =  D/pU2a .  (19) 

where  L  and  D  are  the  lift  and  drag  forces  acting  on  the  cylinder,  respectively,  and  p 
is  the  density  of  the  fluid.  Integrating  the  pressure  and  shear  stresses  over  the  surface 
of  the  cylinder,  we  can  express  these  in  (r,  0)-coordinates  as 

CL  -  CLp+CLf  =  ^£*[-(^  +  <u6]cos0d0  (20) 


The  wake  behind  a  circular  cylinder 


473 


Figure  14.  Equivorticity  contours  for  Re  =  200,  «  =  2.07.  Dashed  (solid)  lines  represent  constant 
positive  (negative)  vorticity  values  with  a  constant  increment  of  0.5,  with  the  magnitude  of  the 
weakest  contour  level  shown  being  0.5.  (a)  t  =  5.0,  (b)  9.0,  (c)  13.0,  ( d )  17.0,  (e)  21.0,  (/)  24.0. 
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Figure  15.  Trajectory  of  the  core  of  the  first  vortex  for  Re  =  200,  a  =  2.07 
obtained  from  equivorticity  contours. 


Figure  16.  Instantaneous  streamlines  for  Re  =  200,  at  =  3.25  at  various  times,  viewed  from  the 
reference  frame  described  in  figure  3.  The  plotting  convention  is  as  in  figure  3.  The  values  of  ^  plotted 
areO,  ±0.1,  ±0.2,  ±0.3,  ±0.4,  ±0.5,  ±0.6,  ±0.7,  ±0.8,  ±0.9,  ±1.0,  with  an  increment  of  ±0.2 
thereafter,  (a)  t  =  5.0,  (b)  9.0,  (c)  16.0,  (d)  24.0. 


Cn  —  CD-  +  Cr 


where  the  subscripts  p  and  /  denote  contributions  due  to  the  pressure  and  friction, 
respectively,  and  the  subscript  b  denotes  quantities  evaluated  on  the  cylinder.  We  note 
that  (20)  and  (21)  differ  from  the  expressions  given  by  Badr  et  al.  (1989)  by  a  sign,  due 
to  a  difference  in  the  definition  of  vorticity.  Figure  21  (a,b)  shows  the  temporal 
evolution  of  the  lift  and  drag  coefficients  at  various  a  for  /  ^  24.  Negative  values  of  CL 


Figure  17.  Comparison  of  computed  (left)  and  experimental  (right)  instantaneous  streamlines  for 
Re  =  200,  a  =  3.25.  Reference  frame  and  camera  motion  are  as  described  in  figure  12.  (a)  t  =  2.0,  (b) 
<.0,  (c)  6.0,  ( d )  8.0. 
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Figure  18.  Equivorticity  contours  for  Re  =  200,  a  =  3.25.  Dashed  (solid)  lines  represent  constant 
positive  (negative)  vorticity  values  with  a  constant  increment  of  0.5.  with  the  magnitude  of  the 
weakest  contour  level  shown  being  0.5.  (a)  t  =  8.0,  ( b )  12.0,  (c)  24.0,  (d)  32.0,  (e)  35.0,  (/)  41.0,  (g) 
48.0,  (h)  54.0. 
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Figure  19.  Instantaneous  streamlines  for  Re  =  200,  a  =  3.25:  (a)  t  =  32.0,  viewed  from  a  frame 
translating  with  the  core  of  the  second  vortex;  (ft)  t  =  54.0,  viewed  from  a  frame  translating  with  the 
core  of  the  third  vortex.  Dashed  (solid)  lines  represent  constant  non-negative  (negative) 
streamfunction  values  with  increments  of  b\/r  =  0.2,  including  f-0.  Note  that  in  (b),  eddy  E,  has 
passed  from  the  field  of  view.  4 


correspond  to  a  lift  force  in  the  negative  ^-direction.  The  time-periodic  nature  of  CL 
is  well  established  for  a  —  0.5  and  1.0.  At  higher  values  of  a,  however,  more  time  is 
required  for  periodicity  to  be  established,  since  the  second  and  subsequent  eddies  form 
and  are  shed  much  later,  as  discussed  above.  In  figure  22 (a,  b),  the  pressure  and  skin 
friction  contributions  to  the  lift  and  drag  are  shown  separately  for  a  =  1 .0.  Similarly 
small  viscous  contributions  to  CL  and  CD  are  found  for  all  other  values  of  a 
investigated.  It  is  clear  that  lift  and  drag  are  largely  due  to  the  pressure  force,  consistent 
with  previous  work  showing  that  the  Magnus  effect  is  primarily  an  inviscid 
phenomenon. 


6.  Discussion 

Our  computations  of  the  temporal  development  of  the  flow  generated  by  a  circular 
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Figure  20.  Trajectory  of  the  core  of  the  first  vortex  for  Re  =  200.  a  =  3.25  obtained  from 

equivorticity  contours. 


cylinder  started  impulsively  from  rest  into  steady  rotatory  and  rectilinear  motion  at 
Re  —  200  show  that,  for  the  largest  value  of  a  (3.25)  investigated  by  Coutanceau  & 
Menard  (1985),  vortex  shedding  continues  after  the  first  vortex  is  shed,  contrary  to 
earlier  conclusions.  It  is  likely  that  these  authors  were  led  to  conclude  that  vortex 
shedding  at  Re  —  200  is  suppressed  for  a  2  because  their  experimental  facility  did  not 
allow  for  visualization  of  the  flow  for  a  long  enough  time.  Even  if  that  limitation  had 
been  overcome,  however,  it  is  likely  that  their  flow  visualization  technique  (which 
approximately  yields  instantaneous  streamlines)  would  have  failed  to  reveal  the 
presence  of  the  second  vortex,  due  to  masking  by  the  large  velocities  induced  in  the 
near  wake  by  the  rapidly  rotating  cylinder. 

We  note  that  for  a  =  3.25  the  time  interval  between  the  shedding  of  the  first  and 
second  vortices  is  much  longer  than  the  time  required  to  shed  the  first.  We  conjecture 
that  for  subsequent  vortices,  the  interval  between  shedding  of  the  2nth  and  (2/i+  l)th 
vortices  will  be  considerably  shorter  than  the  interval  between  shedding  of  the 
(2n+  l)th  and  (2w  +  2)th. 

We  also  observe  that  for  a  =  2.07  and  3.25  the  second  vortex  (shed  from  the 
downstream-moving  side  of  the  cylinder)  is  much  smaller  than  the  first,  and  conjecture 
that,  in  general,  the  2nth  vortex  shed  will  be  significantly  smaller  than  the  (2 n+  l)th. 
Unlike  the  non-rotating  (a  =  0)  case,  there  is  no  requirement  that  the  vortices  shed 
from  opposite  sides  of  the  cylinder  be  of  equal  magnitude,  or  even  that  consecutive 
vortices  be  shed  alternately  from  opposite  sides  of  the  cylinder.  In  fact,  there  does  not 
appear  to  be  any  reason  why  vortices  cannot  be  shed  from  a  single  (upstream-moving) 
side  of  the  cylinder. 

From  a  more  general  standpoint,  we  can  consider  the  computational  and 
experimental  results  in  the  broader  context  of  three-dimensional  flows.  For  this 
purpose,  it  is  useful  to  characterize  the  asymptotic  stability  of  the  flow  in  terms  of  a 
translational  Reynolds  number  defined  by  Re,  =  2 aU/v  and  a  rotational  Reynolds 
number  defined  by  Rer  =  2alQ/v,  equivalent  to  Re  and  a. Re,  respectively.  Thus, 
uniform  flow  past  a  rotating  circular  cylinder  can  be  considered  in  the  (Ret,  Rer)-plane, 
a  quarter  of  which  is  shown  in  figure  23.  From  the  computational  results  of  Jackson 
(1987)  and  Zebib  (1987)  and  the  experimental  work  cited  therein,  we  know  that 
uniform  steady  two-dimensional  flow  past  a  non-rotating  cylinder  (Rer  =  0)  is  stable 
for  Re,  $  45,  at  which  point  a  supercritical  Hopf  bifurcation  to  a  time-periodic  two- 
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-igure  21 .  Temporal  evolution  of  (a)  the  lift  and  (b)  drag  coefficients  for  Re  =  200  and  various  a  for 
0  *£  t  ^  24; - ,  a  =  0.5; . a  =  1.0; - ,  a  =  2.07; - ,  a  =  3.25. 

dimensional  oscillating  wake  flow  occurs.  On  the  other  hand,  for  a  circular  cylinder 
undergoing  steady  rotation  only  ( Ret  =  0),  Walowit,  Tsao  &  DiPrima  (1964)  have 
shown  that  the  steady  purely  azimuthal  flow  {V  =  eff£2a2/r)  is  linearly  stable  for 
Rer  <  11,  at  which  point  it  becomes  unstable  with  respect  to  a  steady  axisymmetric 
flow  consisting  of  Taylor  vortices. 

We  therefore  conjecture  that  steady  two-dimensional  flows  past  a  rotating  circular 
cylinder  are  stable  with  respect  to  infinitesimally  small  but  otherwise  arbitrary 
disturbances  in  a  region  (I  in  figure  23)  of  the  ( Ret ,  /?er)-quarter-plane  including  the 
origin.  As  the  curve  (shown  schematically  and  referred  to  below  as  the  ‘stability 
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Figure  22.  Contributions  of  pressure  and  skin  friction  to  (a)  lift  and  (b)  drag  for  Re  =  200  and 
a.  m  1.0  for  0  <  t  ^  24; - ,  skin  friction; . .  pressure; - ,  total. 

boundary’)  bounding  Region  I  is  crossed,  we  propose  that  for  small  a  (note  that  a  is 
the  slope  of  a  line  connecting  the  origin  to  a  point  in  the  quarter-plane)  the  ensuing 
motion  is  time-periodic  and  two-dimensional,  while  for  large  a  it  is  steady  and  three- 
dimensional  (axisymmetric  in  the  limit  Ret~*  0).  It  would  thus  appear  that 

(a)  there  exists  at  least  one  point  on  the  stability  boundary  at  which  the  nature  of 
the  bifurcation  from  steady  two-dimensional  flow  changes,  and 

(b)  there  exist  two  regions  adjacent  to  I,  in  which  an  unsteady  two-dimensional  flow 
should  be  realizable  (Region  II)  and  in  which  a  steady  three-dimensional  flow  should 
be  realizable  (Region  III). 

Although  it  is  known  for  Rer  =  0  (Williamson  1988)  and  for  Rer  4=  0  (Badr  el  al. 
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Figure  23  Schematic  division  of  the  (Re,,  /te,)-quarter-plane.  Two-dimensional  steady  solutions 
are  stable  in  region  I.  Regions  11  and  Ill  are  described  in  the  text. 


1990)  that  the  unsteady  two-dimensional  flow  becomes  unstable  with  respect  to 
unsteady  three-dimensional  flows  at  sufficiently  high  Ret,  the  bifurcation  structure 
remains  to  be  determined.  Among  the  questions  to  be  answered  are 

(i)  If  Regions  II  and  III  have  a  common  boundary,  how  does  a  transition  occur  from 
unsteady  two-dimensional  flow  to  steady  three-dimensional  flow? 

(ii)  If  Regions  II  and  III  do  not  have  a  common  boundary,  what  lies  between  them? 

(iii)  How  does  the  bifurcation  from  steady  two-dimensional  flow  change  along  the 
stability  boundary? 

From  the  standpoint  of  controlling  laminar  two-dimensional  vortex  shedding  from 
a  circular  cylinder  by  using  either  ‘active’  control  (a  term  that  in  the  fluid  mechanics 
literature  has  come  to  mean  a  time-periodic,  and  frequently  harmonic  input,  e.g. 
S2(t+  T)  —  Q(t)  or  Q(t)  =  Q0  sin  2nft)  or  feedback  control,  the  present  calculations 
show  that  the  nature  of  the  vorticity  field  as  well  as  the  vortex  shedding  process  can  be 
significantly  altered  by  cylinder  rotation.  To  the  best  of  our  knowledge,  the  only 
investigations  to  dat'-  concerning  the  effect  of  a  time-dependent  rotation  rate  on  the 
shedding  process  are  the  experimental  studies  of  Taneda  (1977,  1978)  and  Tokumaru 
&  Dimotakis  (1991),  and  the  combined  experimental  and  theoretical  studies  of 
Okajima,  Takata  &  Asanuma  (1975),  Mo  (1989),  and  Wu,  Mo  &  Vakili  (1989). 

The  experimental  work  of  Taneda  (1978)  for  30  <  Re  <  300  (including  the  range  in 
which  Hopf  bifurcation  (Jackson  1987;  Zebib  1987)  of  the  steady  two-dimensional 
flow  at  Re  %  45  leads  to  unsteady  asymmetric  two-dimensional  solutions  exhibiting 
vortex  shedding)  shows  that  for  sufficiently  large  values  of  the  amplitude  and  frequency 
of  time-harmonic  rotatory  oscillations,  vortex  shedding  (and  indeed  the  formation  of 
attached  vortices)  can  be  eliminated.  Moreover,  Taneda’s  flow  visualizations  (see  his 
figure  3  d)  indicate  that  under  certain  circumstances,  a  flow  can  be  generated  that  is 
nearly  symmetric  about  6  =  0.  This  suggests  that  the  flow  can  be  driven  to  a  symmetric 
state,  and  provides  reason  to  believe  that  it  can  be  stabilized  (in  the  control-theoretical 
sense;  cf.  Kuo  1975)  about  a  symmetric  state  in  which  no  vortex  shedding  occurs. 


The  authors  wish  to  warmly  thank  Dr  M.  Coutanceau  for  providing  the  experimental 


482  Y.-M.  Chen,  Y.-R.  Ou  and  A.  J.  Pearls tein 

results  which  appear  in  this  paper.  This  work  was  supported  by  a  National  Research 
Council  Resident  Research  Associateship  awarded  to  the  first  author.  Support  was 
received  from  the  National  Aeronautics  and  Space  Administration  under  NASA 
Contract  No.  NAS1-18605  and  the  Air  Force  Office  of  Scientific  Research  under 
AFOSR  Grant  89-0079  while  the  second  author  was  in  residence  at  the  Institute  for 
Computer  Applications  in  Science  and  Engineering  (ICASE),  NASA-Langley  Research 
Center.  The  work  of  the  third  author  was  supported  by  NSF  Grants  MSM-8451157 
and  CTS-9017181  and  AFOSR  Grant  90-0156. 


REFERENCES 

Badr,  H.  M.  &  Dennis,  S.  C.  R.  1985  Time-dependent  viscous  flow  past  an  impulsively  started 
rotating  and  translating  circular  cylinder.  J.  Fluid  Mech.  158,  447-488. 

Badr,  H.  M.,  Dennis,  S.  C.  R.  &  Young,  P.  J.  S.  1989  Steady  and  unsteady  flow  past  a  rotating 
circular  cylinder  at  low  Reynolds  numbers.  Computers  &  Fluids  17,  579-609. 

Badr,  H.,  Coutanceau,  M.,  Dennis,  S.  &  Menard,  C.  1986  On  the  phenomenon  of  vortex 
transposition  of  coalescence  in  separated  flows.  C.  R.  Acad.  Sci.  Paris  302(11),  1127-1130. 

Badr.  H.  M„  Coutanceau,  M„  Dennis,  S.  C.  R.  &  Menard,  C.  1990  Unsteady  flow  past  a  routing 
circular  cylinder  at  Reynolds  numbers  103  and  10*.  J.  Fluid  Mech.  220,  459-484. 

Chang,  C.-C.  &  Chern,  R.-L.  1992  Vortex  shedding  from  an  impulsively  surted  rotating  and 
translating  cylinder.  J.  Fluid  Mech.  235,  265-298. 

Chen,  Y.-M.  1989  Numerical  simulation  of  the  unsteady  two-dimensional  flow  in  a  time-dependent 
doubly-connected  domain.  PhD  dissertation.  University  of  Arizona. 

Chen,  Y.-M.,  Koniges,  A.  E.  &  Anderson,  D.  V.  1989  ILUBCG2-1 1 :  Solution  of  1 1-banded 
nonsymmetric  linear  equation  systems  by  a  preconditioned  biconjugate  gradient  routine. 
Comput.  Phys.  Commun.  55,  359-365. 

Collins,  W.  M.  &  Dennis,  S.  C.  R.  1973  Flow  past  an  impulsively  started  circular  cylinder.  J.  Fluid 
Mech.  60,  105-127. 

Coutanceau,  M.  &  Menard,  C.  1985  Influence  of  roution  on  the  near-wake  development  behind 
an  impulsively  started  circular  cylinder.  J.  Fluid  Mech.  158,  399-446. 

Diaz,  F.,  Gavalda,  J.,  Kawall,  J.  G.,  Keffer,  J.  F.  &  Giralt,  F.  1983  Vortex  shedding  from 
a  spinning  cylinder.  Phys.  Fluids  26,  3454-3460. 

Eaton,  B.  E.  1987  Analysis  of  laminar  vortex  shedding  behind  a  circular  cylinder  by  computer-aided 
flow  visualization.  J.  Fluid  Mech.  180,  117-145. 

Fasel,  H.  1976  Investigation  of  the  stability  of  boundary  layers  by  a  finite-difference  model  of  the 
Navier-Stokes  equations.  J.  Fluid  Mech.  78,  355-383. 

Fornberg,  B.  1980  A  numerical  study  of  steady  viscous  flow  past  a  circular  cylinder.  J.  Fluid  Mech. 
98,  819-855. 

Glauert,  M.  B.  1957a  A  boundary  layer  theorem,  with  applications  to  rotating  cylinders.  J.  Fluid 
Mech.  2,  89-99. 

Glauert,  M.  B.  19576  The  flow  past  a  rapidly  routing  circular  cylinder.  Proc.  R.  Soc.  Lond.  A  242, 
108-115. 

Ingham,  D.  B.  1983  Steady  flow  past  a  routing  cylinder.  Computers  &  Fluids  11,  351-366. 

Ingham,  D.  B.  &  Tang,  T.  1990  A  numerical  investigation  into  the  steady  flow  past  a  routing 
circular  cylinder  at  low  and  intermediate  Reynolds  numbers.  J.  Comput.  Phys.  87,  91-107. 

Jackson,  C.  P.  1987  A  finite-element  study  of  the  onset  of  vortex  shedding  in  flow  past  variously 
shaped  bodies.  J.  Fluid  Mech.  182,  23-45. 

Kimura,  T.,  Tsutahara,  M.  &  Wang,  Z.  1992  Wake  of  a  rotating  circular  cylinder.  AIAA  J.  30, 
555-556. 

Klewicki,  J.  C.  &  Falco,  R.  E.  1991  On  accurately  measuring  sutistics  associated  with  small-scale 
turbulent  boundary  layers  using  hot-wire  probes.  J.  Fluid  Mech.  219,  119-142. 

Koromilas,  C.  A.  &  Telionis,  D.  P.  1980  Unsteady  laminar  separation:  an  experimental  study. 
J.  Fluid  Mech.  97,  347-384. 


The  wake  behind  a  circular  cylinder  483 

Krahn,  E.  1955  The  laminar  boundary  layer  on1  a  rotating  cylinder  in  crossflow.  NAVORD  Rep. 
4022. 

Kuo.  B.  1975  Automatic  Control  Systems .  3rd  Edn.  Prentice-Hall. 

Lugt,  H.  J.  1979  The  dilemma  of  defining  a  vortex.  In  Recent  Developments  in  Theoretical  and 
Experimental  Fluid  Mechanics  (ed.  U.  Muller,  K.  G.  Roesner  &  B.  Schmidt),  pp.  309-321. 
Springer. 

Lugt,  H.  J.  &  Haussling,  H.  J.  1974  Laminar  flow  past  an  abruptly  accelerated  elliptic  cylinder  at 
45°  incidence.  /.  Fluid  Mech.  65.  71 1-734. 

Lyulka.  V.  A.  1977  Numerical  solution  of  the  problem  of  the  rotation  of  a  cylinder  in  a  flow  of  a 
viscous  incompressible  fluid.  Vyschisl.  Mat.  mat.  Fiz.  17, 470-480.  (Translated  in  USSR  Comput. 
Maths  and  Math.  Phys.  17(2),  178-188,  1978.) 

Mo,  J.  1989  An  investigation  of  wake  flow  around  a  cylinder  with  rotational  oscillations.  PhD 
dissertation.  University  of  Tennessee  Space  Institute. 

Modi.  V.  J.,  Mokhtarian,  F.  &  Yokomizo,  T.  1990  Effect  of  moving  surfaces  on  the  airfoil 
boundary-layer  control.  J.  Aircraft  27,  42-50. 

Moore.  D.  W.  1957  The  flow  past  a  rapidly  rotating  cylinder  in  a  uniform  stream.  J.  Fluid  Mech. 
2.541-550. 

O Kajima,  A.,  Takata,  H.  &  Asanuma.  T.  1975  Viscous  flow  around  a  rotationally  oscillating 
circular  cylinder.  Rep.  532.  Institute  of  Space  and  Aeronautical  Science,  University  of  Tokyo. 

Orszag,  S.  A.  1980  Spectral  methods  for  problems  in  complex  geometries.  J.  Comput.  Phys.  37, 
70-92. 

Payne,  R.  B.  1958  Calculations  of  unsteady  flow  past  a  circular  cylinder.  J.  Fluid  Mech.  4, 
81-86. 

Perry,  A.  E.,  Chong.  M.  S.  &  Lim.  T.  T.  1982  The  vortex-shedding  process  behind  two-dimensional 
bluff  bodies.  J.  Fluid  Mech.  116.  77-90. 

Prandtl,  L.  1925  The  Magnus  effect  and  windpowered  ships.  Naturwissenschaften  13,  93-108. 

Prandtl,  L.  &  Tietjens,  O.  G.  1934  Applied  Hvdro-  and  Aeromechanics  (transl.  J.  P.  Den  Hartog 
1957).  Dover. 

Reddy,  R.  N.  &  Thompson,  J.  F.  1977  Numerical  solution  of  incompressible  Navier-Stokes 
equations  in  the  integro-differential  formulation  using  boundary-fitted  coordinate  systems. 
AIAA  Paper  77-650. 

Shkadova,  V.  P.  1982  Rotating  cylinder  in  a  flowing  viscous  incompressible  fluid.  Akad.  Nauk 
SSSR ,  Izv.  Mekh.  Zhidk.  Gaza  (1).  16-21.  (Translated  in  Fluid  Dyn.  17  (1),  12-16,  1982.) 

Simuni.  L.  M.  1967  Solution  of  certain  problems  in  flow  of  a  viscous  fluid,  associated  with  a  cylinder 
and  a  sphere.  Izv.  Sib.  Otdel.  Akad.  Nauk  SSSR ,  Ser.  Tekh.  Nauk.  (2),  23-28. 

Swanson,  W.  M.  1961  The  Magnus  effect:  a  summary  of  investigations  to  date.  Trans.  ASME  D: 
J.  Basic  Engng  83,  461-470. 

Ta  Phuoc  Loc  1975  Etude  numerique  de  l’ecoulement  d'un  fluide  visqueux  incompressible  autour 
d’un  cylinder  fixe  ou  en  rotation.  Effet  Magnus.  J.  Mec.  14,  109-134. 

Taneda,  S.  1977  Visual  study  of  unsteady  separated  flows  around  bodies.  Prog.  Aero.  Sci.  17, 
287-348. 

Taneda,  S.  1978  Visual  observations  of  the  flow  past  a  circular  cylinder  performing  a  rotatory 
oscillation.  J.  Phys.  Soc.  Japan  45.  1038-1043. 

Taneda,  S.  1980  Visualization  of  unsteady  flow  separation.  In  Flow  Visualization  II  (ed.  W. 
Merzkirch),  pp.  253-257.  Hemisphere. 

Tang,  T.  &  Ingham,  D.  B.  1991  On  steady  flow  past  a  rotating  circular  cylinder  at  Reynolds 
numbers  60  and  100.  Computers  &  Fluids  19,  217-230. 

Taslim.  M.  E.,  Kinney,  R.  B.  &  Paolino.  M.  A.  1984  Analysis  of  two-dimensional  viscous  flow 
over  cylinders  in  unsteady  motion.  AIAA  J.  22.  586-594. 

Telionis,  D.  P.  1981  Unsteady  Viscous  Flows.  Springer. 

Tennant,  J.  S„  Johnson.  W.  S.  &  Krothapalli,  A.  1976  Rotating  cylinder  for  circulation  control 
on  an  airfoil.  J.  Hydronaut.  10,  102-105. 

Ting,  L.  1983  On  the  application  of  the  integral  invariants  and  decay  laws  of  vorticity  distributions. 
J.  Fluid  Mech.  127,  497-506. 


484  Y.-M.  Chen,  Y.-R.  Ou  and  A.  J.  Pearlstein 

Tokumaru,  P.  T.  &  Dimotakis,  P.  E.  1991  Rotary  oscillation  control  of  a  cylinder  wake.  J.  Fluid 
Mech.  224.  77-90. 

Townsend.  P.  1980  A  numerical  simulation  of  Newtonian  and  visco-elastic  flow  past  stationary  and 
rotating  cylinders.  J.  Non-Newtonian  Fluid  Mech.  6.  219-243. 

Vinokur,  M.  1983  On  one-dimensional  stretching  functions  for  finite-difference  calculations. 
J.  Comput.  Phys.  50.  215-234. 

Walowit.  J..  Tsao,  S.  &  DiPrima.  R.  C.  1964  Stability  of  flow  between  arbitrarily  spaced  concentric 
cylindrical  surfaces  including  the  effect  of  a  radial  temperature  gradient.  Trans.  ASME  E:  J. 
Appl.  Mech.  31,  585-593. 

Wambecq,  A.  1978  Rational  Runge-Kutta  methods  for  solving  systems  of  ordinary  differential 
equations.  Computing  20.  333-342. 

Wang,  C.  M.  &  Wu,  J.  C.  1986  Numerical  solution  of  steady  Navier-Stokes  problems  using  integral 
representations.  AlAA  J.  24,  1305-1312. 

Werle,  H.  1984  Hydrodynamic  visualization  of  the  flow  around  a  streamlined  cylinder  with  suction. 
Cousteau-Malavard  turbine  sail  model.  La  Recherche  Aerospatiale  (English  Edition)  4.  29-38. 

Weston.  R.  P.  &  Liu,  C.  H.  1982  Approximate  boundary  condition  procedure  for  the  two- 
dimensional  numerical  solution  of  vortex  wakes.  AlAA  Paper  82-0951. 

Williamson,  C.  H.  K.  1985  Sinusoidal  flow  relative  to  circular  cylinders.  J.  Fluid  Mech.  155, 
141-174. 

Williamson,  C.  H.  K.  1988  The  existence  of  two  stages  in  the  transition  to  three-dimensionality  on 
a  cylinder  wake.  Phys.  Fluids  31,  3165-3168. 

Wood.  W.  W.  1957  Boundary  layers  whose  streamlines  are  closed.  J.  Fluid  Mech.  2,  77-87. 

Wu,  J.,  Mo,  J.  &  Vakili,  A.  1989  On  the  wake  of  a  cylinder  with  rotational  oscillations.  AlAA  Paper 
89-1024. 

Wu,  J.  C.  1975  Velocity  and  extraneous  boundary  conditions  of  viscous  flow  problems.  AlAA  Paper 
75-47. 

Wu,  J.  C.  1976  Numerical  boundary  conditions  for  viscous  flow  problems.  AlAA  J.  14,  1042-1049. 

Wu,  J.  C.  &  Sampath,  S.  1976  A  numerical  study  of  viscous  flow  around  an  airfoil.  AlAA  Paper  76- 
337. 

Wu,  J.  C.  &  Thompson,  J.  F.  1973  Numerical  solutions  of  time-dependent  incompressible  Navier- 
Stokes  equations  using  an  integro-differential  formulation.  Computers  &  Fluids  1,  197-215. 

Zang,  T.  A.,  Wong,  Y.  S.  &  Hussaini,  M.  Y.  1982  Spectral  multigrid  methods  for  elliptic  equations. 
J.  Comput.  Phys.  48,  485-501. 

Zebib.  A.  1987  Stability  of  viscous  flow  past  a  circular  cylinder.  J.  Engng  Maths  21,  155-165. 


Appendix  B 


On  the  Determination  of  Solenoidal  or  Compressible  Velocity 
Fields  from  Measurements  of  Passive  or  Reactive  Scalars 


Arne  J.  Pearlstein  and  Bonnie  N.  Carpenter 
Department  of  Mechanical  and  Industrial  Engineering 
University  of  Illinois  at  Urbana-Champaign 
1206  West  Green  Street 
Urbana,  IL  61801 


Submitted  to 
Physics  of  Fluids 


Abstract 


Recently,  several  techniques  have  been  proposed  for  the  determination  of  two-  or  three-timensionaJ 
velocity  fields  from  measurements  of  one  passive  scalar.  Here,  we  show  that  measurements  of  one 
scalar  and  knowledge  of  the  equation  governing  its  transport  determine  a  velocity  field  only  up  to  an 
additive  vector  field  locally  perpendicular  to  the  gradient  of  the  scalar  field,  but  otherwise  arbitrary.  We 
then  discuss  procedures  for  selecting  a  unique  velocity  field  from  among  the  uncountable  infinity 
consistent  with  the  scalar  transport  data  and  equation.  We  show  if  that  the  Iterative  procedure  of  Dahm  et 
al.  ( Phys .  Fluids  A  4, 2191-2206, 1992)  for  "solution"  of  a  singular  Inear  equation  system  (obtained  using 
only  measurements  of  the  scalar  and  the  equation  governing  its  transport)  does  in  fact  converge,  that 
outcome  is,  at  best,  an  artifact  of  spatial  discretization,  and  that  as  the  resolution  improves,  the  problem 
becomes  increasingly  ill-conditioned,  becoming  ill-posed  in  the  limit  of  perfect  resolution.  We  then 
develop  a  method  for  determining  an  n-dimensional  (n  ■  2  or  3)  solenoidal  (divergence-free)  velocity  field 
from  measurements  of  n  - 1  passive  or  reactive  scalars.  Finally,  we  show  how  the  velocity  field  in  an  n- 
dimensional  compressible  flow  can  be  determined  from  measurements  of  density  and  n  - 1  passive  or 
reactive  scalars. 
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I.  Introduction 

Determination  of  the  velocity  field  in  any  but  the  simplest  flows  is  a  central  problem  in  experimental 
and  theoretical  fluid  mechanics.  From  an  experimental  standpoint,  it  has  been  possible  for  some  time  to 
accurately  determine  the  time  dependence  of  one  or  more  velocity  components  at  one  or  a  small  number 
of  points,  using  hot-wire,  hot-film,  laser-Doppler  (LDV),  or  other  pointwise  measurement  techniques. 

More  recently,  techniques  have  been  developed  for  estimation  of  one  or  more  velocity  components 
at  more  than  one  point,  using  single-point  methods  in  a  scanning  or  flying  mode  (e.g.,  flying  hot-wire1  or 
scanning  LDV2).  In  “scanning”  mode,  data  is  typically  acquired  at  one  point  at  a  time  along  a  curve  x,(t). 
The  data  are  thus  neither  full-field  nor  simultaneous,  and  are  Smited  to  one  or  at  most  two  velocity 
components. 

In  the  last  ten  yea.s,  several  experimental  techniques  for  estimating  the  velocity  field  have  been 
developed  that  take  advantage  of  essentially  Lagrangian  descriptions  of  the  flow.  These  techniques, 
including  particle  image  velodmetry  (PIV)  and  its  variants3*4,  as  well  as  photochemical,  photochromic. 
and  luminescence  methods5-6  in  which  fluid  elements  are  marked  (e.g.,  using  laser-induced 
photochemistry  or  Raman  excitation  of  vibrational  states),  have  in  common  that  they  attempt  to  follow  the 
trajectories  of  either  neutrally  buoyant  small  particles  or  'marked”  fluid  elements  advected  by  the  flow,  in 
order  to  infer  some  or  all  of  the  components  of  the  velocity  field.  For  the  particle-based  techniques, 
spatial  resolution  is  Bmited  by  several  factors,  including  the  size  of  the  particles  that  can  be  distinguished, 
the  number  of  particles  that  can  be  tracked  from  frame  to  frame  without  ambiguity,  the  volume  fraction 
and  size  of  particles  that  can  be  added  without  seriously  altering  the  flow,  and  the  computational  effort 
required  to  process  the  data.  Other  limitations  of  particle-based  techniques  have  recently  been  discussed 
by  Huang,  Fiedler,  and  Wang7-6.  For  the  photochemical  or  luminescent  techniques,  only  a  fraction  of  the 
volume  of  interest  can  be  "tagged,”  because  if  all  of  the  fluid  is  "maiked,”  then  none  can  be  distinguished. 
For  either  particle  or  marked  fluid  techniques,  data  at  any  time  are  acquired  only  at  the  current  locations 
of  the  particles  or  marked  fluid  elements.  We  refer  to  these  velocity  estimation  techniques  as  being 
"multi-point." 
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Concurrently  with  development  of  multi-point  velocimetry  techniques,  there  has  been  a  rapid 
development  of  optical  techniques  for  the  full-field"  determination  of  temperature  (e.g.,  by  means  of 
Schlieren  or  liquid  crystal  techniques9-10)  or  concentration  (dyes  by  absorption,  other  species  by 
fluorescence  or  phosphorescence,  e.g.,  laser-induced  fluorescence11).  Spatial  resolution  of  these  scalar 
field  measurements  is  limited  only  by  the  resolution  of  the  optics.  This  type  of  data  is  routinely  acquired  at 
high  framing  rates  in  planes  by  rapid  scanning  through  a  three-dimensional  volume.  Thus,  in  two 
dimensions  it  is  possible  to  obtain  simultaneous  full-field  scalar  data,  while  in  three  dimensions,  the 
simultaneity  and  full-field  character  of  the  data  acquired  is  limited  only  by  the  scanning  and  framing  rates. 
We  refer  to  this  data  as  being  simultaneous  and  full-field,  subject  to  imitations  on  temporal  resolution 
imposed  by  framing  rate,  and,  in  three  dimensions,  limitations  on  spatial  resolution  in  the  third  dimension 
imposed  by  scanning  or  framing  rates. 

Recently,  techniques  have  been  advanced12-13  for  determining  velocity  fields  from  full-field  (optical) 
measurements  of  advected  scalars.  Such  techniques  have  the  potential  to  determine  fully  three- 
dimensional  velocity  fields  with  excellent  spatial  resolution,  a  key  advantage  over  PIV-type  methods, 
which  are  ultimately  limited  by  particle  size  and  ioadng  constraints.  In  §11  we  show  that  the  velocity  field 
cannot  be  determined  solely  from  measurements  of  a  single  advected  scalar,  and  that  in  the  pubished 
"scalar  imaging  velocimetry"  procedure  of  Dahm  et  al.12,  the  equations  satisfied  by  the  velocity  field  admit 
an  infinity  of  solutions,  each  of  which  is  indeterminate  by  an  additive  arbitrary  vector  field  locally 
perpendicular  to  the  gracfiertt  of  the  measured  scalar  field.  In  §111,  we  critically  assess  three  techniques 
that  have  been  proposed  for  the  selection  of  one  solution  from  this  uncountable  set.  In  §IV  we  develop  a 
method  for  uniquely  determining  an  n-dimensional  (n  ■  2  or  3)  solenoidal  (divergence-free)  velocity  field 
from  measurements  of  n  - 1  passive  or  reactive  scalars,  making  expfidt  use  of  the  solenoidal  character  of 
the  flow.  In  §V  we  show  how  the  procedure  can  be  extended  to  compressible  flows  for  which  the  density 
field  can  be  measured.  A  discussion  follows  in  §VI. 
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It.  NonunlqiMness  of  Velocity  Fields  Determined  Solely  from  Measurement  of 
One  Scalar  Reid 

Recently,  Dahm  et  al.12  have  proposed  that,  for  a  conserved  passive  scalar  £  whose  distribution 
evolves  according  to 

£  +  u.V£  =  DV*£  .  (1) 


simultaneous  full-field  measurements  of  £  suffice  to  determine  the  velocity  field  u.  Their  discussion 
pertains  specifically  to  three-dmensional  flows  (“four-dimensional  vector  velocity  field  measurements*  in 
the  title  and  elsewhere  refers  to  time  as  a  fourth  cfi  mansion).  Furthermore,  Dahm  et  al.  claim  that  their 
technique  does  not  "explcitly  require  V.u  =  0".  From  that  and  the  claim  that  (1)  and  measurements  of  £ 
suffice  to  determine  u,  it  would  be  logical  to  infer  that  their  technique  is  appicable  to  two-  or  three- 
dimensional  compressible  as  well  as  solenoidal  (i.e.,  V.u=0)  velocity  fields. 

In  fact,  the  method  of  Dahm  et  al.,  in  which  u,,  the  component  of  u  parallel  to  the  gradient  of  the 
scalar  field  (V£),  is  determined  from  measurements  of  £  and 


DV2£-8£/at 


(2) 


and  u  itself  is  said  to  be  determined  from 

Uj  =  u  •  ev; 


(3a) 


and 


Vuj|  =  u  •  Ve^  + V  uT.  ev^  ,  (3b) 

is  inherently  incapable  of  uniquely  determining  the  component(s)  of  the  velocity  field  perpendicular  to  U)t. 
To  see  this,  we  note  that 


u<;=u„ev;  +h(x,t)e1  (4) 

is  a  velocity  field  satisfying  (3a),  where  ex  is  any  unit  vector  locally  orthogonal  to  e?;  and  h(x,t)  is  an 
arbitrary  scalar  function  of  position  and  time.  Hence,  the  solutions  of  (3a, b)  are  neither  discrete  nor 
countable.  Furthermore,  we  note  that  (3b)  can  be  rewritten  as 

Vug  =V(u.eV;)  ,  (5) 
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from  which  it  is  apparent  that  (4)  satisfies  (3b)  also.  Hence,  any  velocity  field  determined  by  solving  (3a,b) 
alone  (using  measurements  of  C)  is  necessarily  nonunique,  because  an  arbitrary  vector  field 
perpendicular  to  VC  can  be  added  to  give  a  Afferent  velocity  field  also  satisfying  (3a,b).  This  is  because 
the  only  relationship  between  the  scalar  field  and  the  velocity  field  is  (1 ),  which  allows  no  distinction  to  be 
made  between  velocity  fields  of  the  form  (4),  so  long  as  ex  is  locally  perpendicular  to  the  gradient  of  the 
scalar.  Put  another  way,  the  component(s)  (one  in  two  dimensions  and  two  in  three  dimensions)  of  the 
velocity  field  that  piay(s)  no  part  in  advecting  the  scalar  cannot  be  uniquely  determined  solely  from 
measurements  of  the  scalar.  This  conclusion  is  independent  of  the  relative  dffusivities  of  momentum  and 
the  scalar. 

Thus,  it  is  clear  that  measurements  of  one  passive  scalar  do  not  suffice  to  determine  an  arbitrary 
velocity  field.  In  the  context  of  machine  vision,  where  in  (1)  C  and  u  represent  the  image  irradiance  and 
optical  flow  field,  respectively,  and  D  =  0  ( correspond ng  to  infinite  Schmidt  number  or  Prandtl  number  in 
the  present  context),  the  inherent  indeterminacy  of  a  vector  field  obtained  solely  from  measurements  of 
the  scalar  field  that  it  advects  is  well  known14'15,  and  is  referred  to  as  the  aperture  problem  or  aperture 
ambiguity.  As  discussed  by  Tokumaru  and  Dimotakis13,  one  or  more  other  "elements"  (e.g.,  one  or  more 
addtional  equations,  minimization  of  some  functional,  etc.)  is  required  in  order  to  render  u  unique.  As 
discussed  below,  we  note  that  the  nature  of  the  additional  element(s)  used  to  select  a  unique  u 
determines  how  well  the  selected  u  corresponds  to  the  actual  velocity  field. 

Alternatively,  the  claim12  that  (3b)  has  a  unique  solution  for  u  because  (3b)  constitutes  three  linear 
equations  in  three  unknowns  is  easily  seen  to  be  false  since  the  right-hand  side  of  (3b)  is  identical  to  the 
right-hand  side  of  (5),  from  which  it  is  easily  seen  that  (4)  satisfies  (5)  (and  hence  3b)  for  any  h(x,t)  so 
long  as  the  vector  field  ex  is  locally  perpendicular  to  VC . 

III.  Selecting  a  Unique  Velocity  Reid 

The  real  issue  then,  as  recognized  by  Tokumaru  and  Dimotakis13,  is  how  one  goes  about  selecting 
the  correct  u  from  among  the  uncountable  infinity  consistent  with  the  scalar  transport  data  and  equation. 
Clearly,  any  technique  that  sufficiently  constrains  the  velocity  field  will  select  a  unique  u.  Thus,  rendering 
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u  unique  is  not  equivalent  to  selecting  the  correct  u.  In  the  following  subsections,  we  discuss  three 
procedures  that  accomplsh,  or  claim  to  accomplsh,  the  selection  of  a  unique  u. 

IB JL  The  Iterative  Inversion"  Approach 

Dahm  et  al.  have  claimed12  that  "solution"  of  (3b)  by  a  particular  iterative  technique  selects  the 
correct  u  (i.e.,  the  velocity  field  that  satisfies  conservation  ot  mass)  from  the  infinity  of  other  solutions  of 
(3a)  and  (3b).  Specifically,  in  reference  12,  (3b)  is  discretized  and  written  as 

b  =  Ai  v(m)+  A*  (6) 

at  each  time,  where  each  element  of  the  vector  b  is  an  approximation  to  a  component  on  the  left-hand 
side  of  (3b)  at  a  point  for  which  scalar  data  are  available,  the  two  terms  on  the  right-hand  side  are 
discretizations  of  the  corresponding  terms  on  the  right-hand  side  of  (3b).  and  v  is  a  vector  of  unknown 
velocity  components  at  the  discrete  points.  It  was  originally  claimed12  that  "Successive  iterations  based 
on  Eq.  (7)"  can  be  made  until  the  velocity  field  converges  to  a  self-consistent  result".  In  addttion  the 
authors  claim  that  "Stable,  rapid,  and  accurate  convergence  is  achieved".  The  correctness  of  the 
computed  velocity  field  was  attributed  to  the  smallness  of  the  second  term  in  (6),  from  which  it  was  said  to 
follow  that  the  initial  iterate  (v*0^  A^b)  was  necessarily  close  to  the  correct  solution  of  (3a, b). 

In  fact,  from  the  indeterminacy  of  the  solutions  of  the  continuous  problem  (3a,b),  it  follows  that  for  a 
faithful  discretization 

b  =Av  (7) 

of  (3b),  the  matrix  A  -A^-t-A^  is  either  singular  (exactly  or  numerically;  the  latter  situation  gives  an  III- 

corxfitioned  linear  equation  system,  so  the  two  possibilities  are  in  practice  equivalent)  at  any  spatial 

resolution,  or  is  singular  in  the  limit  of  infinite  resolution,  in  which  case  the  method  necessarily  breaks 
down  as  the  spatial  resolution  of  the  data  improves  indefinitely.  Here  we  show  that  any  claim  that  the 
iteration  (6)  "converges"  is  inconsistent  with  the  singularity  of  A.  For  general  A,  it  can  be  verified  by 
direct  substitution  that  the  solution  of 


*  equivalent  to  our  (6) 
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b-A1v{m)+(A-A1)v(,n-1) 


(8) 


Since 

XBj=(l-B)-1  (10) 

if  the  infinite  sum  exists,  we  see  that  det(AA~1)  *  0  is  a  necessary  condition  for  the  existence  of  a 
"converged"  formal  solution 

k-1  r  ..j 

v-ttn  Aj’j  [l-AA-1]  b  =  A-,D-(l-AA-1)r1b*A-1(AA-1r1b  .  (11) 

J»o 

of  (8).  But  A  is  singular,  so  no  "converged"  solution  of  (8)  exists,  no  matter  how  A  is  decomposed. 
Thus,  the  meanings  of  the  convergence  claims  made  by  Dahm  et  al.12  are  not  clear.  With  regard  to  the 
statement12  that  "The  ultimate  level  of  improvement  attainable  through  further  iterations  is.  of  course, 
fimited  by  numerical  processing  errors  introduced  by  the  discrete  implementation  of  Eq.  (7)’,  we  note  that 
in  the  limit  of  infinite  resolution  and  no  noise,  the  matrix  A  is  exactly  singular,  and  the  method  necessarily 
fails.  If  in  fact  "converged"  results  are  produced  at  all,  that  would  simply  be  an  artifact  of  the 
discretization.  Moreover,  since  no  solution  of  the  continuous  problem  can  be  obtained  by  the  iterative 
method  described,  any  "solution"  of  the  discrete  problem  obtained  by  this  technique  should  be  expected 
to  depend  sensitively  on  both  noise  in  the  data  and  the  details  of  the  discretization. 


III.B.  The  "Variational”  Approach 

More  recently,  it  has  been  claimed  that  in  three  dimensions  u  can  be  determined  from  measurement 
of  a  single  scalar  if  one  minimizes  the  positive  semidefinite  integral 


.VC-DV2cJ+$(v.u)2  +X(VU;Vu)JJdV  , 


(12) 


where  the  value  of  the  exponent  j  is  variously  given  as  1  or  2,  and  4>  and  X  are  constant  weighting 
factors.  This  is  an  example  of  how  an  additional  "element"  can  be  used  to  select  a  unique  u. 
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It  is  clear  that  for  X  ■  0,  minimization  of  (12)  is  mathematically  equivalent  to  satisfying  (1)  and 
V»u*0.  In  the  following  section,  it  is  shown  that  this  approach  is  correct  for  two-dimensional  sotonoktal 
flows.  For  three-dimensional  flows,  (1)  and  V.u  =  0  are  insufficient  to  uniquely  determine  u.  On  the 
other  hand,  for  X  *0,  one  must  hope  that  minimization  of  (12)  results  in  both  satisfaction  of  (1)  and 
V»u*0,  as  well  as  in  minimization  of  the  viscous  dissipation.  It  has  been  claimed  that  the  viscous 
dissipation  achieves  a  minimum  in  physically  realzed  flows,  and  that  minimization  of  E  selects  the  correct 
u. 

In  fact,  the  class  of  flows  for  which  it  can  be  shown  that  the  viscous  dissipation  is  minimized 
(compared  to  all  other  kinematically  admissible  flows  satisfying  the  same  boundary  conditions)  is  small 
and  rather  uninteresting.  It  consists  of  those  flows  for  which  either 

a)  the  flow  is  steady  and  the  inertial  terms  vanish  identically  (Helmholtz16),  or 

b)  the  viscous  terms  are  derivable  from  a  single-valued  scaiar  potentiaJ  (Rayleigh17). 

Furthermore,  for  sufficiently  large  X,  it  is  likely  that  minimization  of  E  win  drive  the  dissipation  to  zero  at 
the  expense  of  satisfying  (1)  and  V  »u *  0 . 

It  is  clear  that  minimization  of  (12)  cannot  be  expected  to  yield  a  correct  u  for  compressible  flow. 

Neither  condition  a)  nor  b)  obtains  in  the  flows  to  which  the  method  has  been  appfied. 


III.C.  The  "Smoothness"  Approach 

Tokumaru  and  Dimotakis13  have  recently  developed  an  essentially  Lagrangian  technique  that 


considers  situations  for  which 


(13) 


where/  is  the  imaging  resolution,  and  tj-tg  is  the  time  interval  between  consecutive  images  of  the 
scalar  field.  They  integrate  (1)  (with  the  diffusive  term  set  to  zero)  to  get 

CK(x,t1)]-CK(x,t0)]  =  0  ,  (14) 

where  the  displacement  field  is  denoted  by  ,  and  the  initial  locations  of  the  fluid  elements  are  given 


by  $(x,t0)  =  x. 
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Tokumaru  and  Dimotakis  render  the  velocity  field  unique  by  minimizing  the  integral  (over  the 
measurement  volume)  of  the  sum  of  the  squares  of  the  left-hand  side  of  (14)  from  zero,  and  a  measure  of 
smoothness  of  the  displacement  field  £(x,t)  and  its  spatial  derivatives.  This  method  requires  solution  of 
a  global  linear  equation  system  at  each  time  in  order  to  find  a  unique  displacement  field  that  minimizes 
the  error  in  (14)  while  maintaining  its  own  smoothness.  Since  there  is  no  physical  requirement  that  the 
“smoothest*  velocity  field  be  the  correct  one  in  any  particular  flow,  one  cannot  a  priori  assume  that  this 
technique  yields  *he  correct  u.  However,  as  pointed  out  by  Hildreth18  in  the  machine  vision  context, 
although  the  "velocity  field  of  least  variation  is,  in  general,  not  the  physically  correct  one",  it  is  often 
qualitatively  similar.  Tokumaru  and  Dimotakis13  have  appGed  their  method  to  experimental  data  from 
several  carefully  chosen  flows  in  which  they  have  created  and  imaged  appropriate  scalar  distributions. 
Their  results  for  steady  circular  Couette  flow  (with  a  time -dependent  scalar  field)  are  excellent,  while  those 
obtained  using  single  two-dimensional  slices  of  the  three-dimensional  transitional  flow  past  a  circular 
cylinder  are  qualitatively  reasonable. 

IV.  Unique  Determination  of  n-dimensional  Solenoidal  Velocity  Fields  from 
Measurements  of  n-1  Scalar  Fields 

In  this  section,  we  describe  a  technique  employing  full-field  measurements  of  n  - 1  scalars  to 
determine  solenoidal  (divergence-free,  or  "incompressible")  n-dimensional  velocity  fields.  Our  objective  is 
to  develop  a  method  that  is  "consistent",  by  which  we  mean  two  things.  First,  the  method  should  extract 
the  correct  velocity  field  when  there  is  sufficient  scalar  data  to  do  so.  Second,  the  error  in  the  extracted 
velocity  field  should  go  to  zero  as  the  spatial  and  temporal  resolution  of  the  scalar  measurements 
improve.  In  general,  these  requirements  dictate  that  the  additional  "element"  used  to  render  u  unique  be 
satisfied  exactly  in  the  flow,  and  not  implicit  in  the  scalar  transport  equation.  For  two-dmensional  flows, 
we  show  that  measurement  of  one  scalar  (temperature  or  concentration)  suffices  to  uniquely  determine  a 
solenoidal  velocity  field,  whereas  for  three-dimensional  flows,  measurements  of  two  scalars  (temperature 
and  concentration  or  two  independent  concentrations)  having  nonparallel  gradients  are  required.  Our 
technique  is  in  principle  applicable  to  steady  or  unsteady,  laminar  or  turbulent,  two-  or  three-dimensional 
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flows.  Key  features  are  that  no  velocity  components  need  be  measured,  the  number  of  measured  scalars 
required  is  one  less  than  the  spatial  timensionaEty,.  and  the  method  should  give  a  velocity  field  that 
faithfully  approximates  the  correct  u  as  the  spatial  or  temporal  resolution  of  the  scalar  measurements 
improves.  The  computations  are  local  in  a  sense  made  more  precise  below,  and  do  not  require  solution 
of  a  large  linear  (or  nonlinear)  equation  system.  For  steady  or  time-periodic  three-dimensional  flows,  the 
two  scalars  need  not  be  measured  simultaneously. 

Our  technique  is  based  on  the  recognition  that  if  the  transport  equation(s)  for  the  scalar(s)  is  (are) 
linear  in  the  velocity  components,  then  it  (they)  and  the  so-called  incompressible  continuity  equation 

V»u  =  0,  (15) 

constitute  two  (three)  linear  equations  in  two  (three)  velocity  components. 

Our  analysis  will  be  presented  in  a  Cartesian  coordinate  system.  (Results  for  general  curvilinear 
coordinate  systems  mil  be  presented  in  §IV  for  the  compressible  case.)  In  two  dimensions,  we  present 
the  development  in  terms  of  an  equation  for  transport  of  a  single  species  (the  concentration  of  which  is 
denoted  by  C)  in  a  binary  fluid.  We  allow  for  that  species  to  react,  with  kinetics  that  are  essentially 
arbitrary  so  long  as  they  are  known  and  the  rate  of  reaction  depends  solely  on  C.  (In  the  development  of 
the  three-dimensional  analysis,  we  show  how  temperature-dependent  kinetics  can  be  accommodated.) 
Under  a  broad  range  of  conditions,  the  scalar  equation  for  transport  of  C  can  be  approximated  by 

^+u.VC  =  DV2C+R(C),  (16) 

ot 

where  D  is  the  binary  diffusion  coefficient  and  V2  is  the  two-dimensional  Cartesian  Laplacian.  If 
temperature  is  the  passive  scalar  for  which  full-field  measurements  are  available,  then  (16)  will  be 
replaced  by  an  energy  equation,  identical  to  (16)  if  D  is  replaced  by  k  and  the  reaction  rate  R(C)  is  set  to 
zero.  In  two  dimensions,  (16)  can  be  written  as 

+  u  — —  +  v  — —  =  DV2C  +  R(C)  ,  (17) 

3t  3x  dy 

which  together  with  the  two-dimensional  form 
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Equations  (19a)  and  (19b)  are  two  Inear,  first-order,  uncoupled  equations  for  u  and  v. 

In  three  dimensions  we  consider  two  scalars,  at  least  one  of  which  is  a  concentration  variable.  Hie 
second  scalar  can  be  either  temperature  or  another  concentration.  In  what  follows,  we  will  take  the 
second  scalar  to  be  temperature,  and  assume  that  the  thermal  properties  are  temperature-independent 
and  that  the  viscous  dissipation  heating  term  in  the  energy  equation  can  be  neglected.  Only  trivial 
changes  in  the  discussion  are  required  if  both  scalars  are  concentrations.  In  (1 6).  V2  is  now  the  full  three- 
dimensional  Lapladan.  The  two  transport  equations  are  then 


and 


— —  +  u— —  +  v—— •  +  w— — ■  =  DV2C  +  R(T,C), 
dt  3x  3y  3z  v 


£+u£+v|I 

dt  dx  dy 


kV2 


T  + 


(-AH) 

pcp 


R(T,C), 


(21a) 


(2lb) 


which  together  with  the  three-dimensional  version 


of  (15)  constitute  three  linear  equations  for  u,  v,  and  w.  We  note  that  any  kinetic  scheme  (including  one 
with  temperature-dependent  rates)  can  be  accommodated,  so  long  as  the  rate  expression  is  a  known 
expression  of  the  form  R(T.C). 


Three  linear,  first-order  equations 

(VCxVT).V(M/a)=[VTxV(H2/o)-VCxV(H1/a)]  .  e* 

(VC  x  VT) .  V(v/y)  = [VT  x  V^/y)  -  VC  x  V(H,/y)]  .  Cy 

(VC  x  VT) .  V(w/8)  =[VT  x  VfHz/S)  -  VC  x  V(H,/8)]  .  e* 


analogous  to  (19a)  and  (19b)  can  be  derived  from  (21a,b)  and  (22),  where 

H1=DVzC-^+R(T,C) , 


kV2T  -  £ + R(T,C), 

dt  pCp 

and 

=  acar_ac£T 
“  ay  dz~az  ay  * 


acaracar 

dz  dx  dx  dz’ 


and 


acai_acai 

ax  3y  ay  ax 


(23a) 

(23b) 

(23c) 

(24a) 

(24b) 

(25a) 

(25b) 

(25c) 


are  the  components  of  VC  x  VT. 

Equations  (19a,b)  or  (23a-c)  constitute  systems  of  linear  first-order  uncoupled  hyperbolic  equations 
in  the  velocity  components  u  and  v,  and  u,  v,  and  w,  respectively.  Given  the  scalar  field  data,  the  solution 
of  these  equation  systems  is  made  unique  by  specification  of  appropriate  boundary  conditions.  (As  the 
unknowns  are  components  of  the  velocity,  the  required  boundary  conditions  will  involve  velocity 
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components,  rather  than  the  transported  scalars.)  Equations  (19a, b)  or  (23a -c)  can  then  be  integrated 
along  characteristics.  We  note  that  in  two  dimensions,  the  characteristics  are  the  same  for  (19a,b),  and 
that  in  three  dimensions,  the  characteristics  are  the  same  for  (23a-c).  We  note  that  no  time  derivatives  of 
the  velocity  components  appear  in  (19a, b)  or  (23a-c),  so  that  no  initial  data  on  the  velocity  field  are 
required.  A  desirable  consequence  of  the  fact  that  (19a,b)  and  (23a-c)  are  not  initial  value  problems  is 
that  errors  in  the  determination  of  the  velocity  field  do  not  grow  temporally. 

In  two  dimensions,  it  is  dear  that  at  any  point  at  which  either  component  of  the  gradient  of  the  scalar 
vanishes,  (19a)  or  (19b)  will  have  a  singularity.  In  three  dimensions,  equations  (23a),  (23b),  or  (23c)  will 
be  singular  at  points  where  a.  y,  or  S  vanish.  From  the  standpoint  of  extracting  velocity  fields  from 
experimental  scalar  data,  and  more  importantly,  in  order  to  wisely  select  boundary  conditions,  ditfusivity 
ratios,  and  other  factors  that  affect  the  scalar  field(s)  from  which  the  velocity  field  is  to  be  extracted,  it  is 
important  to  understand  the  nature  of  the  singularities,  and  to  develop  techniques  for  dealing  with  them,  if 
posstoie. 

First,  we  classify  the  singularities  of  (19a,b)  and  (23a-c).  We  first  consider  the  two-dimensional 
case,  if  exactly  one  component  of  the  scalar  gradient  vanishes,  then  one  component  of  u  can  be 
determined  from  (19a)  or  (19b),  and  the  other  can  be  found  from  (18).  For  example,  if  dC/dx  =  0,  then 
(19a)  can  be  used  to  find  u  and  (17)  can  be  used  to  determine  v  from 


v = 


Dvzc-ac/dt 


(26) 


In  practice,  it  is  advantageous  to  integrate  (19a)  and  compute  v  from  (26),  and  to  compare  the  values  of  u 
and  v  at  each  point  along  the  common  characteristic  to  the  values  obtained  by  integrating  (19b)  and  using 


pv2c-ac/at 


Thus,  when  one  gets  dose  to  a  curve  on  which  dC/dx  =  0  for  example,  the  results  of  the  above  two 


procedures  will  begin  to  differ,  and  one  can  use  the  results  obtained  from  (19b)  and  (26).  That  the  scalar 
field  and  solenoidal  constraint  still  provide  enough  information  to  determine  u  in  this  case  is  easily  seen, 
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since  the  general  case  can  be  rewritten  in  terms  of  local  coordinates  parallel  and  perpendicular  to  VC. 
characterized  by  unit  vectors  and  ex  ,  such  that  the  component  of  VC  in  the  direction  of  ex  is 
zero. 

When  VC  vanishes,  it  might  seem  from  (17)  that  the  scalar  field  provides  no  information  about  u. 
This  is  indeed  true  if  VC  vanishes  identically  in  some  area,  in  which  case  there  is  insufficient  information 
to  determine  u  in  the  interior  of  that  area.  If,  on  the  other  hand,  VC  vanishes  either  at  an  isolated  point 
(corresponding  to  the  intersection  of  the  curves  along  which  dC/dx  =  0  and  dC/dy  =  0)  or  along  a  curve, 
then  we  can  in  fact  determine  u.  In  two  dimensions,  we  begin  by  using  index  notation  to  rewrite  (17)  as 


dC  „ 
Uj-p-  =  6 
JdXj 


i  =  1.2  . 


(28) 


Taking  the  gradient  of  both  sides  of  (28),  we  obtain 

3*C 


u 


dG 


J  dXjXj  dXj 


1=1,2  . 


(29) 


when  the  gracfient  of  C  vanishes.  Thus,  we  have  two  (ineariy  independent)  nonhomogeneous  eqi  ns 
for  the  components  of  u.  We  find  that 


UsS£ix±iSsl 

CxxCyy— Cjjy 


V  -  Gyc«~G«cxy 
CxxCyy-C^y 


(30a,b) 


We  note  that  (30a, b)  can  be  obtained  from  (17)  only  if  VC  =  0.  That  the  soienoidal  constraint  is  not 
required  in  order  to  specify  the  velocity  field  in  this  case  is  due  to  the  fact  that  another  constraint,  namely 
VC  =  0,  obtains. 


Thus,  we  see  that  unless  VC  vanishes  identically  over  an  area,  "singularities"  in  (I9a,b)  are 
’removable*  in  the  sense  that  the  velocity  field  can  be  uniquely  determined  at  points  or  along  curves 
where  one  or  both  components  of  VC  vanish. 

In  three  dimensions,  VC  x  VT  plays  the  part  of  VC,  and  the  vanishing  of  one  or  two  components  of 
the  cross  product  of  the  grarfients  of  the  scalars  corresponds  to  a  change  to  a  local  coortfnate  system 
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with  unit  vectors  parallel  and  tangential  to  the  cross  product.  In  that  case,  formulas  analogous  to  (26)  and 
(27)  can  be  derived  for  the  components  of  u.  If  all  three  components  of  VC  x  VT  vanish  (i.e.,  the 
gradients  of  C  and  T  are  parallel)  in  a  volume,  then  the  scalar  transport  equations  and  the  soienokfal 
constraint  provide  insufficient  information  for  determination  of  u  in  that  volume,  in  the  event  that 
VC  x  VT  vanishes  on  a  cun/e  or  on  a  surface,  expressions  for  u,  v,  and  w  analogous  to  (30a, b)  for  the 
two-dimensional  case  can  be  derived.  The  apparent  singularity  in  (23a-c)  is  completely  removable  unless 
VCxVT  vanishes  in  a  volume. 

Finally,  we  note  that  one  can  easily  show  that  characteristics  originating  on  a  boundary  will  emanate 
into  the  flow  (as  opposed  to  remaining  on  the  boundary)  if  the  normal  component  of  the  scalar  gradient  (in 
two  dimensions)  or  the  normal  component  of  the  cross  product  ol  the  scalar  gradients  (in  three 
dimensions)  is  nonzero.  In  two  dimensions,  this  point  argues  strongly  for  use  of  temperature  as  the 
scalar,  since  most  walls  are  impermeable  to  most  species  (i.e.,  n»  VC  ■  0).  In  three  dimensions,  a 
nonisothermal  wall  through  which  the  second  scalar  does  not  penetrate  will  be  consistent  with 
n  •  (VC  x  VT)  *  0 . 

V.  Unique  Determination  of  n-dimensional  Compressible  Velocity  Fields  from 
Measurements  of  Density  and  n-1  Other  Scalar  Fields 

In  this  section,  we  show  how  the  approach  of  §111  can  be  extended  to  deal  with  the  case  in  which  the 
equation  for  conservation  of  mass  is 

^£  +  u.Vp  +  pV.u  =  0  .  (31) 

dt 

We  again  use 

~  +  u • VC  =  DVZC  +  R(T,C)  (32a) 

3t 

and 

pcJ  — +  u.VT|  =  V.kVT  +  ^^-[^  +  u.Vp|  +  (-AH)R(T,C)  ,  (32b) 

Vdt  ;  P  V*  ) 
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where  0  is  the  thermal  expansion  coefficient,  to  eliminate  two  components  of  u  before  substituting  into 
the  continuity  equation  (this  time  (31)]  to  get  a  hypeiboic  equation  for  the  third  component.  Doing  this 
three  times,  we  get,  in  general  curviinear  coo rtS nates, 


(33) 


where  L=  VC x(VT+  rVp),  r  =  (1-Cp/cv)/(p p),  hi,  ha,  and  h3  are  metric  coefficients,  f- hi hah3, 
Hi  is  as  defined  in  §111,  and 


Ha= 


V.kVT 

PCp 


(-AH) 

P°P 


R(T,C)  . 


(34) 


In  the  compressible  case,  singularity  corresponds  to  the  vanishing  of  components  of  L,  and  the 
approach  is  similar  to  that  described  in  §111. 

For  a  constant  density  fluid,  Cp= cv ,  and  we  recover  the  solenoidal  results  given  in  §111. 

VI.  Discussion 

Determination  of  an  n-dmensional  velocity  field  (n  *  2  or  3)  requires  n  "elements,’  each  constituting 
an  independent  constraint  on  the  velocity  components.  In  our  approach,  for  a  two-dimensional  solenoidal 
velocity  field  these  are  equations  (17)  and  (18),  where  the  derivatives  of  the  scalar  field  in  (17)  are 
obtained  from  experiment.  In  three  dimensions,  we  take  the  three  elements  to  be  equations  (2la,b)  and 
(22),  where  the  scalar  fields  and  their  derivatives  are  obtained  from  experiment.  For  compressible  flow, 
one  can  use  the  continuity  equation  (31)  and  n  - 1  of  equations  (32a)  and  (20b),  where  the  density  and 
advected  scalar(s)  are  obtained  from  experiment. 

We  note  that,  for  solenoidal  u,  one  modification  of  the  approach  presented  above  might  be  to 
employ  the  same  elements  to  compute  one  velocity  component  using  the  corresponding  hyperbolic 
equation,  and  to  use  the  n  - 1  transport  equations  (17,  or  21  a, b),  which  are  linear  and  algebraic  in  the 
other  components,  to  determine  the  remainder  of  the  velocity  field. 
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A  rather  different  set  of  n  elements  might  be  full-field  measurements  of  n-r  velocity  components, 
r  - 1  transport  equations,  and  a  continuity  equation. 

In  any  case,  the  fact  that  n  independent  equations  are  required  to  determine  n  unknowns  means  that 
a  single  transport  equation  along  with  measurement  of  the  corresponding  scalar  field  does  not  uniquely 
determine  a  solenoidal  velocity  field  in  three  dimensions,  or  a  compressible  velocity  field  in  two 
dimensions,  let  alone  a  compressible  velocity  field  in  three  dimensions. 

We  note  that  the  feasibilty  of  our  technique  will  depend  in  large  part  on  one's  ability  to  create  scalar 
distributions  for  which  VC  (in  two  dimensions)  or  VC  x  VT  (in  three  dimensions)  is  not  only  not  zero 
everywhere  in  the  interior  of  an  area  (in  two  dimensions)  or  a  volume  (in  three  dimensions),  but  is  also  not 
small  (compared  to  the  measurement  sensitivity  divided  by  the  spatial  resolution)  in  the  interior  of  a 
volume  or  area.  This  will,  in  general,  require  some  care  to  be  exercised  in  setting  up  the  scalar  field(s). 
In  addition  to  the  choice  of  boundary  conditions  for  the  scalar,  and  some  Rmited  flexibiBty  in  the  choice  of 
diffusivities,  one  may  also  want  to  take  advantage  of  the  fact  that,  since  the  technique  is  applicable  to 
scalars  whose  transport  is  described  by  (16)  in  two  dimensions  and  (19a,b)  in  three  dimensions,  one  can 
consider  various  "volumetric  methods  for  influencing  the  scalar  fields,  in  which  R(C)  or  R(T,C)  is  replaced 
or  augmented  by  a  spatially  nonuniform  "source  term”. 

Finally,  we  note  that  our  technique  is  an  ideal  cantidate  for  parallelization.  The  integration  along 
characteristics  uses  only  local  information  along  the  characteristic,  and  does  not  involve  solution  of  a 
"global"  equation  system,  as  is  required  by  the  "iterative  inversion"  and  "variational"  approaches. 
Moreover,  integration  along  the  characteristics  is  a  process  that  can  be  parceled  out  on  the  basis  of  one 
processor  per  characteristic.  The  fact  that  no  large  linear  equation  system  need  be  solved  (the 
"variational  technique"  described  in  §II.B.  involves  a  linear  equation  system  so  large  as  to  require  iterative 
solution),  makes  the  process  eminently  suitable  for  real-time  use. 
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